ALGO: implement the quad tree

Algorithms/RangeSearch/Core.hs

@@ -0,0 +1,191 @@
+module Algorithms.RangeSearch.Core
+    (quadTree,
+    quadTreeSquares,
+    qtFoldl,
+    qtFoldr,
+    goQuad,
+    findNeighbor,
+    lookupByPath',
+    rootNode,
+    testArr,
+    Orient,
+    Quad,
+    QuadTree)
+  where
+
+import Algebra.VectorTypes
+import Algebra.Vector
+import Data.Foldable (foldlM)
+import Data.Maybe (fromJust)
+import Diagrams.TwoD.Types
+
+
+-- |The quad tree structure.
+data QuadTree a
+  -- |An empty node.
+  = TNil
+  -- |A leaf containing some value.
+  | TLeaf a
+  -- |A node with four children.
+  | TNode (QuadTree a) (QuadTree a)   -- NW NE
+          (QuadTree a) (QuadTree a)   -- SW SE
+  deriving (Show, Eq)
+
+-- |Represents a Quadrant in the 2D plane.
+data Quad = NW | NE
+          | SW | SE
+
+-- |A Crumb used for the QuadTree Zipper.
+data Crumb a = NWCrumb (QuadTree a) (QuadTree a) (QuadTree a)
+             | NECrumb (QuadTree a) (QuadTree a) (QuadTree a)
+             | SWCrumb (QuadTree a) (QuadTree a) (QuadTree a)
+             | SECrumb (QuadTree a) (QuadTree a) (QuadTree a)
+  deriving (Show, Eq)
+
+-- |A list of Crumbs.
+type Breadbrumbs a = [Crumb a]
+
+-- |Zipper for the QuadTree.
+type Zipper a = (QuadTree a, Breadbrumbs a)
+
+-- |Orientation.
+data Orient = North | South | East | West
+
+
+-- |Get a sub-square of the current square, e.g. nw, ne, sw or se.
+nwSq, neSq, swSq, seSq :: Square -> Square
+nwSq ((xl, xu), (yl, yu)) = (,)  (xl, (xl + xu) / 2)  ((yl + yu) / 2, yu)
+neSq ((xl, xu), (yl, yu)) = (,)  ((xl + xu) / 2, xu)  ((yl + yu) / 2, yu)
+swSq ((xl, xu), (yl, yu)) = (,)  (xl, (xl + xu) / 2)  (yl, (yl + yu) / 2)
+seSq ((xl, xu), (yl, yu)) = (,)  ((xl + xu) / 2, xu)  (yl, (yl + yu) / 2)
+
+
+-- |Check whether the current Node is an nw, ne, sw or se child of it's
+-- parent.
+isNWchild, isNEchild, isSWchild, isSEchild :: Zipper a -> Bool
+isNWchild (_, NWCrumb {}:_) = True
+isNWchild _                 = False
+isNEchild (_, NECrumb {}:_) = True
+isNEchild _                 = False
+isSWchild (_, SWCrumb {}:_) = True
+isSWchild _                 = False
+isSEchild (_, SECrumb {}:_) = True
+isSEchild _                 = False
+
+
+-- |Builds a quadtree of a list of points which recursively divides up 2D
+-- space into quadrants, so that every leaf-quadrant stores either zero or one
+-- point.
+quadTree :: [PT]        -- ^ the points to divide
+         -> Square      -- ^ the initial square around the points
+         -> QuadTree PT -- ^ the quad tree
+quadTree pts' sq' = go (flip filter pts' . inRange $ sq') sq'
+  where
+    go []   _ = TNil
+    go [pt] _ = TLeaf pt
+    go pts sq = TNode (quadTree pts . nwSq $ sq) (quadTree pts . neSq $ sq)
+                      (quadTree pts . swSq $ sq) (quadTree pts . seSq $ sq)
+
+
+-- |Get all squares of a quad tree.
+quadTreeSquares :: Square      -- ^ the initial square around the points
+                -> QuadTree PT -- ^ the quad tree
+                -> [Square]    -- ^ all squares of the quad tree
+quadTreeSquares sq (TNil)              = [sq]
+quadTreeSquares sq (TLeaf _)           = [sq]
+quadTreeSquares sq (TNode nw ne sw se) =
+  quadTreeSquares (nwSq sq) nw ++ quadTreeSquares (neSq sq) ne ++
+  quadTreeSquares (swSq sq) sw ++ quadTreeSquares (seSq sq) se
+
+
+-- |Left fold over the tree leafs.
+qtFoldl :: (a -> b -> a) -> a -> QuadTree b -> a
+qtFoldl _ sv (TNil)              = sv
+qtFoldl f sv (TLeaf a)           = f sv a
+qtFoldl f sv (TNode nw ne sw se) = foldl (qtFoldl f) sv [nw, ne, sw, se]
+
+
+-- |Right fold over the tree leafs.
+qtFoldr :: (b -> a -> a) -> a -> QuadTree b -> a
+qtFoldr f sv qt = qtFoldl (\g b x -> g (f b x)) id qt sv
+
+
+-- |Go to nw, ne, sw or se from the current node, one level deeper.
+goNW, goNE, goSW, goSE :: Zipper a -> Maybe (Zipper a)
+goNW (TNode nw ne sw se, bs) = Just (nw, NWCrumb ne sw se:bs)
+goNW _                       = Nothing
+goNE (TNode nw ne sw se, bs) = Just (ne, NECrumb nw sw se:bs)
+goNE _                       = Nothing
+goSW (TNode nw ne sw se, bs) = Just (sw, SWCrumb nw ne se:bs)
+goSW _                       = Nothing
+goSE (TNode nw ne sw se, bs) = Just (se, SECrumb nw ne sw:bs)
+goSE _                       = Nothing
+
+
+-- |Go to the given Quad from the current Node, one level deeper.
+goQuad :: Quad -> Zipper a -> Maybe (Zipper a)
+goQuad q = case q of
+  NW -> goNW
+  NE -> goNE
+  SW -> goSW
+  SE -> goSE
+
+
+-- |Go up to the parent node, if any.
+goUp :: Zipper a -> Maybe (Zipper a)
+goUp (qt, NWCrumb ne sw se:bs) = Just (TNode qt ne sw se, bs)
+goUp (qt, NECrumb nw sw se:bs) = Just (TNode nw qt sw se, bs)
+goUp (qt, SWCrumb nw ne se:bs) = Just (TNode nw ne qt se, bs)
+goUp (qt, SECrumb nw ne sw:bs) = Just (TNode nw ne sw qt, bs)
+goUp _                         = Nothing
+
+
+-- |Get the root node.
+rootNode :: Zipper a -> Zipper a
+rootNode (qt, []) = (qt, [])
+rootNode z        = rootNode . fromJust . goUp $ z
+
+
+-- |Look up a node by a given path of Quads.
+lookupByPath' :: [Quad] -> QuadTree a -> Maybe (Zipper a)
+lookupByPath' qs qt = foldlM (flip goQuad) (qt, []) qs
+
+
+-- |Find the north, south, east or west neighbor of a given node.
+findNeighbor :: Orient -> Zipper a -> Maybe (Zipper a)
+findNeighbor ot zr = case ot of
+  North -> go isSWchild isSEchild isNWchild goNW goNE goSW goSE zr
+  South -> go isNWchild isNEchild isSWchild goSW goSE goNW goNE zr
+  East  -> go isNWchild isSWchild isNEchild goNE goSE goNW goSW zr
+  West  -> go isNEchild isSEchild isNWchild goNW goSW goNE goSE zr
+  where
+    go _ _ _ _ _ _ _  (_, []) = Nothing
+    go is1 is2 is3 go1 go2 go3 go4 z@(_, _:_)
+      | is1 z = goUp z >>= go1
+      | is2 z = goUp z >>= go2
+      | otherwise = checkParent
+                      . go is1 is2 is3 go1 go2 go3 go4
+                      . fromJust
+                      . goUp
+                      $ z
+      where
+        checkParent (Just (z'@(TNode {}, _)))
+          | is3 z     = go3 z'
+          | otherwise = go4 z'
+        checkParent (Just z') = Just z'
+        checkParent _         = Nothing
+
+
+
+testArr :: [PT]
+testArr = [p2 (200.0, 450.0),
+           p2 (400.0, 350.0),
+           p2 (100.0, 300.0),
+           p2 (25.0 , 350.0),
+           p2 (225.0, 225.0),
+           p2 (400.0, 150.0),
+           p2 (300.0, 100.0),
+           p2 (300.0, 300.0),
+           p2 (300.0, 350.0),
+           p2 (50.0 , 450.0),
+           p2 (100.0, 25.0)]

Algorithms/RangeSearch/UB4_Punkte.obj

@@ -0,0 +1,11 @@
+v 200.0 450.0
+v 400.0 350.0
+v 100.0 300.0
+v 225.0 225.0
+v 400.0 150.0
+v 300.0 100.0
+v 300.0 300.0
+v 300.0 350.0
+v 100.0 25.0
+v 25.0 350.0
+v 50.0 450.0

Algorithms/RangeSearch/UB4_Suchbereich.obj

@@ -0,0 +1,5 @@
+v 50.0 100.0
+v 250.0 100.0
+v 250.0 475.0
+v 50.0 475.0
+f 1 2 3 4

Algorithms/RangeSearch/UB4_alles_zusammen.obj

@@ -0,0 +1,14 @@
+
+v 50.0 100.0
+v 250.0 100.0
+v 250.0 475.0
+v 50.0 475.0
+f 1 2 3 4
+v 200.0 450.0
+v 400.0 350.0
+v 100.0 300.0
+v 225.0 225.0
+v 400.0 150.0
+v 300.0 100.0
+v 100.0 25.0
+

CG2.cabal

@@ -54,7 +54,7 @@ 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 QueueEx Algorithms.PolygonIntersection.Core
+  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 Algorithms.PolygonIntersection.Core Algorithms.RangeSearch.Core
 
   -- LANGUAGE extensions used by modules in this package.
   -- other-extensions:

GUI/gtk2.glade

@@ -767,7 +767,8 @@ Public License instead of this License.
             <property name="items" translatable="yes">Show points
 Show convex hull
 Show polygons
-Show polygons intersection</property>
+Show polygons intersection
+Show quad tree squares</property>
           </widget>
           <packing>
             <property name="expand">False</property>

Graphics/Diagram/Gtk.hs

@@ -18,6 +18,10 @@ diag p obj@(Object _)
   | alg p == 1 =
       mkDiag (mconcat [convexHPText, convexHP, convexHLs, coordPoints, plotterBG])
              p obj
+  | alg p == 4 =
+      mkDiag (mconcat [squares, coordPointsText, coordPoints, polyLines,
+                       plotterBG])
+             p obj
   | otherwise = mempty
 diag p objs@(Objects _)
   | alg p == 2 =

Graphics/Diagram/Plotter.hs

@@ -4,11 +4,13 @@ module Graphics.Diagram.Plotter where
 
 import Algebra.VectorTypes
 import Algorithms.ConvexHull.GrahamScan
+import Algorithms.RangeSearch.Core
 import Algorithms.PolygonIntersection.Core
 import Data.Monoid
 import Diagrams.Backend.Cairo
 import Diagrams.Prelude hiding ((<>))
 import Graphics.Diagram.Types
+import Graphics.Gloss.Data.Extent
 
 
 -- |Creates a Diagram that shows the coordinates from the points
@@ -33,6 +35,7 @@ pointToTextCoord :: PT -> Diagram Cairo R2
 pointToTextCoord pt =
   text ("(" ++ (show . trim') x ++ ", " ++ (show . trim') y ++ ")") # scale 10
   where
+    trim' :: Double -> Double
     trim' x' = (fromInteger . round $ x' * (10^2)) / (10.0^^2)
     (x, y) = unp2 pt
 
@@ -165,6 +168,26 @@ convexHStepsLs = GifDiag chs
           (strokeTrail . fromVertices $ vt') # moveTo (head vt') # lc col
 
 
+-- FIXME: hardcoded dimensions
+squares :: Diag
+squares = Diag f
+  where
+    f p (Object []) = mempty
+    f p (Object vt) =
+      mconcat
+      $ (\((xmin, xmax), (ymin, ymax)) -> rect (xmax - xmin) (ymax - ymin) #
+        moveTo (p2 (
+                     ((xmax + xmin) / 2),
+                     ((ymax + ymin) / 2)
+                   )
+               ) #
+        lw ultraThin)
+      <$> (quadTreeSquares ((0,500), (0,500)) . quadTree vtf $ (dX p, dY p))
+      where
+        vtf = filterValidPT p vt
+    f _ _ = mempty
+
+
 -- |Creates a Diagram that shows an XAxis which is bound
 -- by the dimensions given in xD from DiagProp.
 xAxis :: Diag