POLYINT: first try of polygon intersection algorithm in O(n)

Algebra/Vector.hs

@@ -3,7 +3,10 @@
 module Algebra.Vector where
 
 import Algebra.VectorTypes
+import Control.Applicative
 import Diagrams.TwoD.Types
+import Graphics.Gloss.Geometry.Line
+import GHC.Float
 import MyPrelude
 
 
@@ -72,6 +75,27 @@ det a b c =
     (cx, cy) = unp2 c
 
 
+-- |Computes the determinant of 2 points.
+det' :: PT -> PT -> Double
+det' a b =
+  ax * by - ay * bx
+  where
+    (ax, ay) = unp2 a
+    (bx, by) = unp2 b
+
+
+-- |Get the point where two lines intesect, if any.
+intersectSeg' :: Segment -> Segment -> Maybe PT
+intersectSeg' (a, b) (c, d) =
+  glossToPt <$> intersectSegSeg (ptToGloss a)
+                                (ptToGloss b)
+                                (ptToGloss c)
+                                (ptToGloss d)
+  where
+    ptToGloss = (\(x, y) -> (double2Float x, double2Float y)) <$> unp2
+    glossToPt = p2 . (\(x, y) -> (float2Double x, float2Double y))
+
+
 -- |Get the orientation of 3 points which can either be
 -- * clock-wise
 -- * counter-clock-wise
@@ -95,3 +119,26 @@ notcw a b c = case getOrient a b c of
 -- |Sort X and Y coordinates lexicographically.
 sortedXY :: [PT] -> [PT]
 sortedXY = fmap p2 . sortLex . fmap unp2
+
+
+-- |Apply a function on the coordinates of a point.
+onPT :: ((Double, Double) -> (Double, Double)) -> PT -> PT
+onPT f = p2 . f . unp2
+
+
+-- |Compare the y-coordinate of two points.
+ptCmpY :: PT -> PT -> Ordering
+ptCmpY p1' p2' = compare ((snd . unp2) p1') ((snd . unp2) p2')
+
+
+-- |Compare the y-coordinate of two points.
+ptCmpX :: PT -> PT -> Ordering
+ptCmpX p1' p2' = compare ((fst . unp2) p1') ((fst . unp2) p2')
+
+
+posInfPT :: PT
+posInfPT = p2 (read "Infinity", read "Infinity")
+
+
+negInfPT :: PT
+negInfPT = p2 (negate . read $ "Infinity", negate . read $ "Infinity")

Algebra/VectorTypes.hs

@@ -5,11 +5,13 @@ module Algebra.VectorTypes where
 import Diagrams.TwoD.Types
 
 
-type Vec = R2
-type PT  = P2
-type Coord  = (Double, Double)
+type Vec     = R2
+type PT      = P2
+type Coord   = (Double, Double)
+type Segment = (PT, PT)
 
 
 data Alignment = CW
                | CCW
                | CL
+  deriving (Eq)

Algorithms/PolygonIntersection/Core.hs

@@ -0,0 +1,180 @@
+module Algorithms.PolygonIntersection.Core where
+
+
+import           Algebra.Vector
+import           Algebra.VectorTypes
+import           Data.Dequeue (BankersDequeue)
+import qualified Data.Dequeue           as Q
+import           Data.List
+import           Data.Maybe
+import           Diagrams.TwoD.Types
+import           MyPrelude
+import           QueueEx
+import           Safe
+
+
+-- |Describes a point on the convex hull of the polygon.
+-- In addition to the point itself, both it's predecessor and
+-- successor are saved for convenience.
+data PolyPT =
+  PolyA {
+    id'   :: PT
+    , pre :: PT
+    , suc :: PT
+    }
+  | PolyB {
+    id'   :: PT
+    , pre :: PT
+    , suc :: PT
+    }
+  deriving (Show, Eq)
+
+
+-- |Shift a list of sorted convex hull points of a polygon so that
+-- the first element in the list is the one with the highest y-coordinate.
+-- This is done in O(n).
+sortLexPoly :: [PT] -> [PT]
+sortLexPoly ps = maybe [] (`shiftM` ps) (elemIndex (yMax ps) ps)
+  where
+    yMax       = foldl1 (\x y -> if ptCmpY x y == GT then x else y)
+
+
+-- | Sort the points of two polygons according to their y-coordinates,
+-- while saving the origin of that point. This is done in O(n).
+sortLexPolys :: ([PT], [PT]) -> [PolyPT]
+sortLexPolys (pA'@(_:_), pB'@(_:_)) =
+  queueToList . go (Q.fromList . sortLexPoly $ pA') $
+                   (Q.fromList . sortLexPoly $ pB')
+  where
+    -- Start recursive algorithm, each polygon is represented by a Queue.
+    -- Traverse predecessor and successor and insert them in the right
+    -- order into the resulting queue.
+    -- We start at the max y-coordinates of both polygons.
+    go :: BankersDequeue PT -> BankersDequeue PT -> BankersDequeue PolyPT
+    go pA pB
+      -- Nothing to sort.
+      | Q.null pA && Q.null pB
+        = Q.empty
+      -- Current point of polygon A is higher on the y-axis than the
+      -- current point of polygon B, so insert it into the resulting
+      -- queue and traverse the rest.
+      -- remark: we don't handle y1 = y2
+      | ptCmpY (fromMaybe negInfPT . Q.first $ pA)
+               (fromMaybe posInfPT . Q.first $ pB) == GT
+        = Q.pushFront
+            (go (maybeShift . snd . Q.popFront $ pA) pB)
+            (PolyA (fromJust . Q.first $ pA)
+              (pre' pA' pA)
+              (suc' pA' pA))
+
+      -- Same as above, except that the current point of polygon B
+      -- is higher.
+      | otherwise
+        = Q.pushFront
+            (go pA (maybeShift . snd . Q.popFront $ pB))
+            (PolyB (fromJust . Q.first $ pB)
+              (pre' pB' pB)
+              (suc' pB' pB))
+
+    pre' xs = fromJust . polySuccessor xs . uQfirst
+    suc' xs = fromJust . polyPredecessor xs . uQfirst
+
+    -- Compare the first and the last element of the queue according
+    -- to their y-coordinate and shift the queue (if necessary) so that
+    -- the element with the highest value is at the front.
+    maybeShift :: BankersDequeue PT -> BankersDequeue PT
+    -- remark: we don't handle y1 = y2
+    maybeShift q = if ptCmpY (fromMaybe posInfPT . Q.first $ q)
+                             (fromMaybe negInfPT . Q.last  $ q) == GT
+                     then q
+                     else shiftQueueRight q
+sortLexPolys _ = []
+
+
+-- |Get the successor of a point on a convex hull of a polygon.
+-- Returns Nothing if the point is not on the convex hull. This
+-- is done in O(n).
+polySuccessor :: [PT] -> PT -> Maybe PT
+polySuccessor pts pt = case index of
+  Nothing      -> Nothing
+  Just index'  -> if index' == (length pts - 1)
+                    then pts `atMay` 0
+                    else pts `atMay` (index' + 1)
+  where
+    index = elemIndex pt pts
+
+
+-- |Get the predecessor of a point on a convex hull of a polygon.
+-- Returns Nothing if the point is not on the convex hull. This
+-- is done in O(n).
+polyPredecessor :: [PT] -> PT -> Maybe PT
+polyPredecessor pts pt = case index of
+  Nothing     -> Nothing
+  Just index' -> if index' == 0
+                   then pts `atMay` (length pts - 1)
+                   else pts `atMay` (index' - 1)
+  where
+    index = elemIndex pt pts
+
+
+-- |Get all points that intersect between both polygons. This is done
+-- in O(n).
+intersectionPoints :: [PolyPT] -> [PT]
+intersectionPoints []  = []
+intersectionPoints xs' =
+  rmdups
+    . (++) (segIntersections . scanLine $ xs')
+    $ intersectionPoints (tail xs')
+  where
+    -- Get the scan line or in other words the
+    -- Segment pairs we are going to check for intersection.
+    scanLine :: [PolyPT] -> ([Segment], [Segment])
+    scanLine xs = (segmentsA xs, sgementsB xs)
+
+    -- Gets the actual intersections between the segments of
+    -- both polygons we currently examine. This is done in O(1)
+    -- since we have max 4 segments.
+    segIntersections :: ([Segment], [Segment]) -> [PT]
+    segIntersections (a@(_:_), b@(_:_))
+      = catMaybes
+          . fmap (\[x, y] -> intersectSeg' x y)
+          $ combinations a b
+    segIntersections _ = []
+
+    -- Gets all unique(!) combinations of two arrays. Both arrays
+    -- are max 2, so this is actually O(1) for this algorithm.
+    combinations :: [a] -> [a] -> [[a]]
+    combinations xs ys = concat . fmap (\y -> fmap (\x -> [y, x]) xs) $ ys
+
+    segmentsA :: [PolyPT] -> [Segment]
+    segmentsA sp@(_:_) = case a of
+      Nothing -> []
+      Just x  -> [(id' x, suc x), (id' x, pre x)]
+      where
+        a = listToMaybe . filter (\x -> case x of
+                                           PolyA {} -> True
+                                           _       -> False) $ sp
+    segmentsA _ = []
+
+    sgementsB :: [PolyPT] -> [Segment]
+    sgementsB sp@(_:_) = case b of
+      Nothing -> []
+      Just x  -> [(id' x, suc x), (id' x, pre x)]
+      where
+        b = listToMaybe . filter (\x -> case x of
+                                           PolyB {} -> True
+                                           _       -> False) $ sp
+    sgementsB _ = []
+
+
+testArr :: ([PT], [PT])
+testArr = ([p2 (200.0, 500.0),
+            p2 (0.0, 200.0),
+            p2 (200.0, 100.0),
+            p2 (400.0, 300.0)],
+
+           [p2 (350.0, 450.0),
+            p2 (275.0, 225.0),
+            p2 (350.0, 50.0),
+            p2 (500.0, 0.0),
+            p2 (450.0, 400.0)])

Algorithms/PolygonIntersection/UB2_einfachCCW.obj

@@ -0,0 +1,15 @@
+v 200.0 500.0
+v 0.0 200.0
+v 200.0 100.0
+v 400.0 300.0
+
+v 350.0 450.0
+v 275.0 225.0
+v 500.0 0.0
+v 350.0 50.0
+v 450.0 400.0
+
+f 1 2 3 4
+f 5 6 7 8
+
+/* ([(200.0, 500.0), (0.0 200.0), (200.0 100.0), (400.0 300.0)], [(350.0 450.0), (275.0 225.0), (500.0 0.0), (450.0 400.0)]) */

Algorithms/PolygonIntersection/UB2_einfachCW.obj

@@ -0,0 +1,10 @@
+v 200.0 500.0
+v 0.0 200.0
+v 200.0 100.0
+v 400.0 300.0
+v 350.0 450.0
+v 275.0 225.0
+v 500.0 0.0
+v 450.0 400.0
+f 4 3 2 1
+f 8 7 6 5 

CG2.cabal

@@ -54,13 +54,13 @@ executable Gtk
   main-is:             GtkMain.hs
 
   -- Modules included in this executable, other than Main.
-  other-modules:       MyPrelude GUI.Gtk Graphics.Diagram.Gtk Graphics.Diagram.Types Graphics.Diagram.Plotter Parser.Meshparser Parser.Core System.FileSystem.FileExt Algebra.Vector Algorithms.ConvexHull.GrahamScan
+  other-modules:       MyPrelude GUI.Gtk Graphics.Diagram.Gtk Graphics.Diagram.Types Graphics.Diagram.Plotter Parser.Meshparser Parser.Core System.FileSystem.FileExt Algebra.Vector Algorithms.ConvexHull.GrahamScan QueueEx
 
   -- LANGUAGE extensions used by modules in this package.
   -- other-extensions:
 
   -- Other library packages from which modules are imported.
-  build-depends:       base >=4.6 && <4.8, diagrams-lib >=1.2 && <1.3, diagrams-cairo >=1.2 && <1.3, transformers >=0.4 && <0.5, glade >=0.12 && <0.13, gtk >=0.12 && <0.13, directory >=1.2 && <1.3
+  build-depends:       base >=4.6 && <4.8, diagrams-lib >=1.2 && <1.3, diagrams-cairo >=1.2 && <1.3, transformers >=0.4 && <0.5, glade >=0.12 && <0.13, gtk >=0.12 && <0.13, directory >=1.2 && <1.3, dequeue >= 0.1.5, multiset-comb >= 0.2.1, gloss >= 1.2.0.1, safe >= 0.3.8
 
   -- Directories containing source files.
   -- hs-source-dirs:
@@ -74,13 +74,13 @@ executable Gif
   main-is:             GifMain.hs
 
   -- Modules included in this executable, other than Main.
-  other-modules:       MyPrelude Graphics.Diagram.Gif Graphics.Diagram.Types Graphics.Diagram.Plotter Parser.Meshparser Parser.Core System.FileSystem.FileExt Algebra.Vector Algorithms.ConvexHull.GrahamScan
+  other-modules:       MyPrelude Graphics.Diagram.Gif Graphics.Diagram.Types Graphics.Diagram.Plotter Parser.Meshparser Parser.Core System.FileSystem.FileExt Algebra.Vector Algorithms.ConvexHull.GrahamScan QueueEx
 
   -- LANGUAGE extensions used by modules in this package.
   -- other-extensions:
 
   -- Other library packages from which modules are imported.
-  build-depends:       base >=4.6 && <4.8, diagrams-lib >=1.2 && <1.3, diagrams-cairo >=1.2 && <1.3, transformers >=0.4 && <0.5, JuicyPixels >= 3.1.7.1
+  build-depends:       base >=4.6 && <4.8, diagrams-lib >=1.2 && <1.3, diagrams-cairo >=1.2 && <1.3, transformers >=0.4 && <0.5, JuicyPixels >= 3.1.7.1, dequeue >= 0.1.5, multiset-comb >= 0.2.1, gloss >= 1.2.0.1, safe >= 0.3.8
 
   -- Directories containing source files.
   -- hs-source-dirs:

MyPrelude.hs

@@ -78,3 +78,10 @@ dupLast xs = xs ++ [last xs]
 if' :: Bool -> a -> a -> a
 if' True  x _ = x
 if' False _ y = y
+
+
+-- |Shift a list k places.
+shiftM :: Int -> [a] -> [a]
+shiftM _ [] = []
+shiftM 0 xs = xs
+shiftM n xs = drop n xs ++ take n xs

QueueEx.hs

@@ -0,0 +1,28 @@
+module QueueEx where
+
+import           Control.Applicative
+import           Data.Dequeue (BankersDequeue)
+import qualified Data.Dequeue           as Q
+import           Data.Maybe
+
+
+-- |Shift a queue to the left, such as:
+-- [1, 2, 3] -> [2, 3, 1]
+shiftQueueLeft :: BankersDequeue a -> BankersDequeue a
+shiftQueueLeft  = (\(b, nq) -> Q.pushBack nq (fromJust b)) <$> Q.popFront
+
+
+-- |Shift a queue to the right, such as:
+-- [1, 2, 3] -> [3, 1, 2]
+shiftQueueRight :: BankersDequeue a -> BankersDequeue a
+shiftQueueRight = (\(b, nq) -> Q.pushFront nq (fromJust b)) <$> Q.popBack
+
+
+-- |Convert a Queue back to a list.
+queueToList :: BankersDequeue a -> [a]
+queueToList q = Q.takeFront (Q.length q) q
+
+
+-- |Unsafe version of Q.first.
+uQfirst :: BankersDequeue a -> a
+uQfirst = fromJust . Q.first