diff --git a/Algebra/Vector.hs b/Algebra/Vector.hs index 8f01ff8..1b5d026 100644 --- a/Algebra/Vector.hs +++ b/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") diff --git a/Algebra/VectorTypes.hs b/Algebra/VectorTypes.hs index a465d91..d6934bc 100644 --- a/Algebra/VectorTypes.hs +++ b/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) diff --git a/Algorithms/PolygonIntersection/Core.hs b/Algorithms/PolygonIntersection/Core.hs new file mode 100644 index 0000000..cc84def --- /dev/null +++ b/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)]) diff --git a/Algorithms/PolygonIntersection/UB2_einfachCCW.obj b/Algorithms/PolygonIntersection/UB2_einfachCCW.obj new file mode 100644 index 0000000..af9770f --- /dev/null +++ b/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)]) */ diff --git a/Algorithms/PolygonIntersection/UB2_einfachCW.obj b/Algorithms/PolygonIntersection/UB2_einfachCW.obj new file mode 100644 index 0000000..283c3c7 --- /dev/null +++ b/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 diff --git a/CG2.cabal b/CG2.cabal index 8d1d28d..e457765 100644 --- a/CG2.cabal +++ b/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: diff --git a/MyPrelude.hs b/MyPrelude.hs index a3a1326..2b05698 100644 --- a/MyPrelude.hs +++ b/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 diff --git a/QueueEx.hs b/QueueEx.hs new file mode 100644 index 0000000..1446204 --- /dev/null +++ b/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