diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,29 @@
+Copyright (c) 2009, Corey O'Connor
+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 the author nor the names of its 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,5 @@
+module Main where
+import Distribution.Simple
+
+main = defaultMain
+
diff --git a/data-spacepart.cabal b/data-spacepart.cabal
new file mode 100644
--- /dev/null
+++ b/data-spacepart.cabal
@@ -0,0 +1,40 @@
+Name:           data-spacepart
+Version:        0.1.1
+License:        BSD3
+License-File:   LICENSE
+Author:         Corey O'Connor <coreyoconnor@gmail.com>
+Maintainer:     Corey O'Connor <coreyoconnor@gmail.com>
+Homepage:       http://www.tothepowerofdisco.com/repo/data-spacepart/
+Package-URL:    http://www.tothepowerofdisco.com/repo/data-spacepart/
+Category:       Data
+Build-Type:     Simple
+Synopsis:       Space partition data structures. Currently only a QuadTree.
+Stability:      alpha
+Description:
+    Space partition data structures. Currently only a QuadTree.
+    .
+    darcs get --partial http:\/\/code.haskell.org\/data-spacepart\/
+    .
+    TODO:
+    .
+        lots.
+    .
+    See README: http:\/\/code.haskell.org\/data-spacepart\/README
+    
+Extra-Source-Files: test/QuadTreeVisualize.hs
+                    test/run_visualize
+                    test/Render.hs
+                    test/run_verify
+                    test/Verify.hs
+                    test/Verify/Data/AABB.hs
+                    test/Verify/Data/QuadTree.hs
+
+Cabal-Version:  >= 1.6
+
+library
+    hs-source-dirs:         src
+    build-depends:          base, vector-space == 0.5.*, mersenne-random == 0.1.*
+    exposed-modules:        Data.QuadTree
+                            Data.AABB
+                            System.Random.Utils
+
diff --git a/src/Data/AABB.hs b/src/Data/AABB.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/AABB.hs
@@ -0,0 +1,193 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE IncoherentInstances #-}
+{-# LANGUAGE UndecidableInstances #-}
+{-# LANGUAGE TypeSynonymInstances #-}
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE NoMonomorphismRestriction #-}
+module Data.AABB
+    ( Boundary(..)
+    , HasBoundary(..)
+    , encloses
+    , intersects
+    )
+    where
+
+import Data.VectorSpace
+import Data.List (foldl')
+
+type Vertex2 a = (a, a)
+vec2 :: Double -> Double -> Vertex2 Double
+vec2 x y = (x, y)
+
+vx :: VectorSpace (v, v) => (v, v) -> v
+vx (x, _) = x
+
+vy :: VectorSpace (v, v) => (v, v) -> v
+vy (_, y) = y
+
+type Edge2 a = (Vertex2 a, Vertex2 a)
+type LineSegment = Edge2 Double
+
+-- "intersects" is a commutative binary predicate on two shapes. 
+class Intersectable s0 s1 where
+   intersects :: s0 -> s1 -> Bool
+
+--instance Intersectable s0 s1 => Intersectable s1 s0 where
+--   intersects s1 s0 = intersects s0 s1
+
+intersections e es = filter (intersects e) es
+
+-- | A 2D axis aligned square.
+-- The boundary_corner defines the lower bound.
+-- The boundary_size is the length of any edge of the square.
+--
+-- The boundary is inclusive on the low extent and exclusive on the max extent.
+--
+-- Used to represent both the 
+-- 0. 2D axis aligned minimum bounding square of an element.
+--
+-- 1. The boundary of a quadtree element
+--
+data Boundary = Boundary
+    {
+        boundary_corner :: Vertex2 Double,
+        boundary_size   :: Double
+    }
+    deriving (Eq, Show)
+
+-- Boundaries b0 and b1 intersect if the min extent of the intersection of b1 with (the plane +x
+-- including b0.p unioned with the plane +y including b0.p) is within b0.
+instance Intersectable Boundary Boundary where
+    intersects b0 b1 = 
+        let c = (MinExtentPlanes $ boundary_corner b0)
+        in if b1 `intersects` c
+            then let (Boundary p _) = intersection b1 c
+                 in intersects b0 p
+            else False
+
+newtype MinExtentPlanes = MinExtentPlanes (Vertex2 Double)
+    deriving (Eq, Show)
+
+-- A boundary intersects the min extent planes if the far extent of the boundary is within the range
+-- defined by the min extent planes.  The comparison is > and not >= since the far extent is the
+-- point just beyond the boundary. Which needs to be just inside the planes in order for the
+-- boundary to be inside the planes.
+instance Intersectable Boundary MinExtentPlanes where
+    intersects b (MinExtentPlanes (min_x, min_y)) =
+        let ((b_min_x, b_min_y), (b_max_x, b_max_y)) = boundary_extents b
+        in if b_min_x == min_x && b_min_y == min_y
+            then True
+            else (b_max_x > min_x) && (b_max_y > min_y)
+
+intersection :: Boundary -> MinExtentPlanes -> Boundary
+intersection (Boundary p size) (MinExtentPlanes min_p) = Boundary (ext_max min_p p) size
+
+instance Intersectable Boundary LineSegment where
+   intersects b l@(p0, p1) =
+-- If any point of the line segment is contained in the boundary then the line segment intersects the
+-- element.
+       intersects b p0 || intersects b p1
+-- If niether point is in the element the line segment could still intersect the boundary. The line
+-- segment must, in this case, intersect an edge of the boundary.
+       || any (intersects l) (boundary_edges b)
+
+--The equations for line intersection are pulled from 
+--  http://local.wasp.uwa.edu.au/~pbourke/geometry/lineline2d/
+-- without much thought.
+
+instance Intersectable LineSegment LineSegment where
+   intersects (p0a, p0b) (p1a, p1b) = 
+       let x1 = vx p0a
+           y1 = vy p0a
+           x2 = vx p0b
+           y2 = vy p0b
+           x3 = vx p1a
+           y3 = vy p1a
+           x4 = vx p1b
+           y4 = vy p1b
+           div = (y4 - y3)*(x2 - x1) - (x4 - x3)*(y2 - y1)
+       in if div < 1e-9 
+           then False
+           else
+           let t0n = (x4 - x3)*(y1 - y3) - (y4 - y3)*(x1 - x3)
+               t0 = t0n / div
+               t1n = (x2 - x1)*(y1 - y3) - (y2 - y1)*(x1 - x3)
+               t1 = t1n / div
+           in t0 > 0.0 && t0 < 1.0 && t1 > 0.0 && t1 < 1.0
+
+union_boundaries :: Boundary -> Boundary -> Boundary
+union_boundaries b0 b1 =
+    let (min0, max0) = boundary_extents b0
+        (min1, max1) = boundary_extents b1
+        p = ext_min min0 min1
+        ext = ext_max max0 max1
+        (w,h) = ext ^-^ p
+        size = max w h
+    in Boundary p size
+
+ext_min (x0,y0) (x1,y1) = (min x0 x1, min y0 y1)
+ext_max (x0,y0) (x1,y1) = (max x0 x1, max y0 y1)
+
+--instance Show Boundary where
+--    show (Boundary p size) = show p ++ " -> " ++ show size
+
+instance Intersectable Boundary (Vertex2 Double) where
+    intersects bounds (px, py) =
+        let (x, y) = boundary_corner bounds
+            s = boundary_size bounds
+        -- If the point is equal to the corner point then consider it intersecting.
+        -- The inclusive nature of the min extent "wins out" over the exclusive nature of the max
+        -- extent.
+        in if x == px && y == py
+            then True
+            else px < (x + s) && px >= x && py < (y + s) && py >= y
+
+{- | A instance of HasBoundary has an axis aligned boundign square defined that entirely encloses
+ - the space represented by the type.
+ -}
+class HasBoundary s where
+    boundary_points :: s -> [Vertex2 Double]
+    boundary_edges :: s -> [Edge2 Double]
+    boundary_edges s = 
+        let ps@(p0 : ps') = boundary_points s
+        in zip ps (ps' ++ [p0])
+    boundary_extents :: s -> (Vertex2 Double, Vertex2 Double)
+    boundary_extents s =
+        let (p0 : ps) = boundary_points s
+            initial_min_extent = p0
+            initial_max_extent = p0
+            union_extents ((min_x, min_y), (max_x,max_y)) (x, y) =
+                let min_x' = min min_x x
+                    min_y' = min min_y y
+                    max_x' = max max_x x
+                    max_y' = max max_y y
+                in ((min_x', min_y'), (max_x', max_y'))
+        in foldl' union_extents (initial_min_extent, initial_max_extent) ps
+    boundary_square :: s -> Boundary
+    boundary_square s =
+        let (min_extent, max_extent) = boundary_extents s
+            width  = fst max_extent - fst min_extent
+            height = snd max_extent - snd min_extent
+            size = max width height
+        in Boundary (fst min_extent, snd min_extent) size
+
+-- A boundary cleary has itself as it's boundary.
+instance HasBoundary Boundary where
+    boundary_points (Boundary p s) = 
+        [ p
+        , p ^+^ (0, s)
+        , p ^+^ (s, s)
+        , p ^+^ (s, 0)
+        ]
+    boundary_extents (Boundary p s) = (p, p ^+^ (s,s))
+    boundary_square b = b
+
+{-| Returns true if the first boundary entirely encloses the second boundary.
+ - This is expected to be reflexive.
+ -}
+encloses :: Boundary -> Boundary -> Bool
+encloses (Boundary (x0,y0) s0) (Boundary (x1,y1) s1) = (x0 <= x1 && x0 + s0 >= x1 + s1) && (y0 <= y1 && y0 + s0 >= y1 + s1)
+
diff --git a/src/Data/QuadTree.hs b/src/Data/QuadTree.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/QuadTree.hs
@@ -0,0 +1,308 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE IncoherentInstances #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE GADTs #-}
+module Data.QuadTree where
+import Data.AABB
+
+import Data.Maybe
+import Data.List ( sortBy )
+import qualified Data.List as List
+import Data.VectorSpace
+
+-- | A 2D binary hierarchical space subdivision of a region. 
+-- All elements contained in the quadtree are required to have a Boundary. This is an axis aligned
+-- box with congruent sides.
+--
+-- Each node of the quadtree is composed of:
+-- 
+-- 0. A list of elements who's shape can be queried for intersection with the quad.  These are all
+-- the elements with a boundary that are fully enclosed by the boundary of this node but not fully
+-- enclosed by a quadrant of this node. 
+--
+-- 1. The Boundary of this node.
+--
+-- 2. The child nodes of this node. Each is a quadrant of this nodes boundary.
+--
+data QuadTree e where
+    QuadTree :: (HasBoundary e)
+                => [e]
+                -> Boundary
+                -> ( Maybe (QuadTree e) , Maybe (QuadTree e)
+                   , Maybe (QuadTree e) , Maybe (QuadTree e)
+                   ) 
+                -> QuadTree e
+
+elements :: QuadTree e -> [e]
+elements (QuadTree es _ _) = es
+
+children :: QuadTree e -> ( Maybe (QuadTree e) , Maybe (QuadTree e)
+                          , Maybe (QuadTree e) , Maybe (QuadTree e)
+                          )
+children (QuadTree _ _ c) = c
+
+instance HasBoundary (QuadTree e) where
+    boundary_points (QuadTree _ bounds _) = boundary_points bounds
+    boundary_edges (QuadTree _ bounds _) = boundary_edges bounds
+    boundary_extents (QuadTree _ bounds _) = boundary_extents bounds
+    boundary_square (QuadTree _ bounds _) = bounds
+
+data Quadrant = 
+        NPQuad | PPQuad
+      | NNQuad | PNQuad
+    deriving (Eq, Show)
+
+{- An element of a quadtree can intersect the boundary of multiple nodes in the quadtree. This
+ - only associates an element with a single node.  This permits the property: forall p. p <: paths
+ - to a leaf node, forall e <: elements in the universe => p will enounter no more than one element
+ - that references e.
+ -
+ - Which, I think, simplifies things. Maybe?
+ - EG:
+ -  let qt = QuadTree.empty
+ -      qt' = QuadTree.insert qt e0
+ -      qt'' = QuadTree.insert q' e1
+ -  In the case where e1 entirely encompasses qt' there would be greater sharing betwee qt'' and qt'
+ -  than if each node in the tree contained references to all elements that intersect that node.
+ -
+ -  On the other hand the query "All elements intersecting this child node of this quadtree." would
+ -  require a full descent from the root to collect the list of elements. I could see this being a
+ -  useful query. 
+ -
+ -  I think this is resolvable. The query necessitates a cursor like structure: The reference to a
+ -  specific child node in a quadtree. Which could transparently cache the parent node element
+ -  references.
+ -} 
+
+pp_quad (QuadTree _ _
+         (  _, mq, 
+            _,  _
+         )
+        ) 
+    = mq
+pn_quad (QuadTree _ _
+         (  _,  _, 
+            _, mq
+         )
+        ) 
+    = mq
+nn_quad (QuadTree _ _
+         (  _,  _, 
+           mq,  _
+         )
+        ) 
+    = mq
+np_quad (QuadTree _ _
+         ( mq,  _, 
+            _,  _
+         )
+        ) 
+    = mq
+
+map_child :: (Maybe (QuadTree e) -> Maybe (QuadTree e)) 
+             -> Quadrant 
+             -> (Maybe (QuadTree e), Maybe (QuadTree e)
+                ,Maybe (QuadTree e), Maybe (QuadTree e)
+                )
+             -> (Maybe (QuadTree e), Maybe (QuadTree e)
+                ,Maybe (QuadTree e), Maybe (QuadTree e)
+                )
+map_child f NPQuad ( np_c, pp_c
+                   , nn_c, pn_c
+                   ) = ( f np_c, pp_c
+                       , nn_c  , pn_c
+                       )
+map_child f PPQuad ( np_c, pp_c
+                   , nn_c, pn_c
+                   ) = ( np_c, f pp_c
+                       , nn_c, pn_c
+                       )
+map_child f NNQuad ( np_c, pp_c
+                   , nn_c, pn_c
+                   ) = ( np_c  , pp_c
+                       , f nn_c, pn_c
+                       )
+map_child f PNQuad ( np_c, pp_c
+                   , nn_c, pn_c
+                   ) = ( np_c, pp_c
+                       , nn_c, f pn_c
+                       )
+
+non_empty_children q =
+    let (np_c, pp_c, nn_c, pn_c) = children q
+    in catMaybes [np_c, pp_c, nn_c, pn_c]
+
+{- | Returns an empty QuadTree without a specific boundary. The default bounds are centered around
+ - (0,0) with a size of 2
+ -
+ - TODO: Alternatively an empty quadtree could have no defined bounds. The bounds would then be
+ - defined on the first insertion. 
+ -}
+empty :: HasBoundary e => QuadTree e
+empty = QuadTree [] (Boundary (-1,-1) 2) empty_children
+
+{- | Returns an empty QuadTree with the given bounds.
+ - The given bounds cannot have a size of 0. This will error out on that case.
+ -
+ - TODO: The user may find it easier for this to accept a 0 sized boundary which is transparently
+ - changed to a non-0 sized boundary on insert.
+ -}
+empty_with_bounds :: HasBoundary e => Boundary -> QuadTree e
+empty_with_bounds (Boundary _ 0.0) = error "Cannot construct a quadtree with 0 sized boundary."
+empty_with_bounds bounds = QuadTree [] bounds empty_children
+
+empty_children = ( Nothing, Nothing
+                 , Nothing, Nothing
+                 ) 
+
+singleton_child NPQuad q = ( Just q , Nothing
+                           , Nothing, Nothing
+                           )
+singleton_child PPQuad q = ( Nothing, Just q
+                           , Nothing, Nothing
+                           )
+singleton_child NNQuad q = ( Nothing, Nothing
+                           , Just q , Nothing
+                           )
+singleton_child PNQuad q = ( Nothing, Nothing
+                           , Nothing, Just q
+                           )
+
+{-| Inserts the given element into the quadtree. 
+ - This inserts the element into a this node or a child quadrant node if the current node encloses
+ - the element.  Otherwise this inserts the element into a new node that is a parent of the given
+ - node.
+ -}
+insert :: (HasBoundary e) => e -> QuadTree e -> QuadTree e
+insert e q =
+    if (boundary_square q) `encloses` (boundary_square e)
+        then insert_self_or_child e q
+        else insert_via_parent e q
+
+{-| Inserts the given element into either a child node of the current node if one of the quadrants
+ - encloses the element.
+ - Otherwise the element is added to the current node's list of elements.
+ -}
+insert_self_or_child :: (HasBoundary e) => e -> QuadTree e -> QuadTree e
+insert_self_or_child e q@(QuadTree es bounds quadrants) =
+    case filter (\(cqb, _) -> cqb `encloses` (boundary_square e)) (quadrant_bounds q) of
+        [child]      -> insert_child child e q 
+        _            -> QuadTree (e : es) bounds quadrants
+
+quadrant_bounds :: QuadTree e -> [(Boundary, Quadrant)]
+quadrant_bounds (QuadTree _ (Boundary p size) _) = 
+    let child_size = size / 2
+        nn_p = p 
+        np_p = p ^+^ (0         , child_size)
+        pp_p = p ^+^ (child_size, child_size)
+        pn_p = p ^+^ (child_size, 0         )
+    in map (\(p, q) -> (Boundary p child_size, q))
+        [ (nn_p, NNQuad)
+        , (np_p, NPQuad)
+        , (pp_p, PPQuad)
+        , (pn_p, PNQuad)
+        ]
+
+{- insert_via_parent adds the given element to a new quadtree, q_e, that is connected to the given
+ - quadtree, q, through a parent tree, q_root. 
+ -
+ - The two quadtrees q and q_e are both children on some path from q_root.
+ -
+ - There is at least one path from q_root to q and q_e. There may be multiple paths?
+ -  let q = (-1, -1) -> 1
+ -      q_e = (0,0) -> 1
+ -      q_root = (-1,-1) -> 2
+ - In the above case there is only one possible q_root with minimum bounds. However there are multiple
+ - mays to connect q and q_e through a parent node.
+ -      q_p_0 = (-2, -2) -> 2 [PP => q]
+ -      q_p_1 = (0,0) -> 2 [NN => q_e]
+ -      q_root = (-2, -2) -> 4 [NN => q_p_0, PP => q_p_1]
+ -
+ - I'm not really sure of how to optimally introduce a node for q_e and connect them through a
+ - parent node.  There are incorrect methods. EG: Always picking the parent quadtree such that the
+ - given quadtree is at a fixed position.  This could result in a search for a new encompasing
+ - parent that never converges.
+ -
+ - The method used here is to add parent nodes to q until a parent node is found that encompass e.
+ - This is a breadth first search of the generated graph
+ - Nodes are parent quadtrees containing q as a child and encompasing e
+ - Edges are directional (q_u, q_v). Each edge represents the operation of adding a parent to q_u
+ - such that q_u is a specific quadrant of the parent.
+ -
+ - Given quadtree q and an element e:
+ -  There is an edge from q for each of PNQuad, PPQuad, NPQuad, NNQuad to a parent quadtree with q
+ -  as the given quadrant.
+ -  This parent quadtree can be generated from q and the quadrant identifier.
+ -}
+
+-- | Adds the element to quadtree via a parent node to the given quadtree.
+-- The parent to add e to is then the first of the possible parents nodes that enclose e.
+insert_via_parent :: (HasBoundary e) 
+                    => e
+                    -> QuadTree e 
+                    -> QuadTree e
+insert_via_parent e q = 
+    let q_root = first (\pq ->  (boundary_square pq) `encloses` (boundary_square e)) (parent_trees q)
+    in insert_self_or_child e q_root
+    where first f = fromJust . List.find f
+
+-- | parent_trees generates all possible parent trees of the given tree (Without memoization) in the
+-- order suitable for a breadth first search.
+parent_trees q = parent_trees' [q]
+    where 
+        parent_trees' (q : qs) = 
+            let parents = imm_parents q
+            in parents ++ parent_trees' (qs ++ parents)
+        imm_parents q_child = map (quadtree_with_child_in_quad q_child) [PNQuad, PPQuad, NPQuad, NNQuad]
+
+quadtree_with_child_in_quad q@(QuadTree _ (Boundary (child_x,child_y) child_size) _) quad 
+    | quad == NPQuad = QuadTree [] (Boundary (child_x, child_y - child_size) parent_size) $ singleton_child quad q
+    | quad == PPQuad = QuadTree [] (Boundary (child_x - child_size, child_y - child_size) parent_size) $ singleton_child quad q
+    | quad == PNQuad = QuadTree [] (Boundary (child_x - child_size, child_y) parent_size) $ singleton_child quad q
+    | quad == NNQuad = QuadTree [] (Boundary (child_x, child_y) parent_size) $ singleton_child quad q
+    where parent_size = child_size * 2
+
+{- I wonder if there is a closed form solution to the search performed by insert_via_parent
+ -
+ -  For all Integer i =>
+ -      The size of the quadrants at this level are equal to 
+ -          size_i = base_size * 2^i
+ -      For all Integer u,v => 
+ -          The corner points of the quadrants are given by
+ -          ( base_point.x + size_i * u, base_point.y + size_i * v)
+ -  The search is for an (i,u,v) such that the quadrant identified by (i,u,v) completely encompases
+ -  the element being inserted.
+ -  For a given i it is possible to find a quadrant that either encompasses the element or
+ -  intersects the elements boundary.
+ -}
+
+-- | Inserts the element in the child identified by the given boundary and Quadrant.
+-- If there is no child at the given quadrant then a child is added and the element is inserted into
+-- the new child.
+insert_child :: (HasBoundary e) 
+                => (Boundary, Quadrant) 
+                -> e 
+                -> QuadTree e
+                -> QuadTree e
+insert_child (cb, quad) e q@(QuadTree es b cs) = 
+    let update_child = Just . insert_self_or_child e . maybe (QuadTree [] cb empty_children) id
+    in QuadTree es b $ map_child update_child quad cs
+
+{- | Returns all elements with boundaries that intersect the given boundary
+ - By case:
+ -  Boundary does not intersect quadtree
+ -  Boundary intersects the quadtree
+ -      All elements at the level of the quadtree could intersect the boundary. Test each element
+ -      for intersection. 
+ -      Descend into the quadrants
+ -}
+query :: (HasBoundary e) => Boundary -> QuadTree e -> [e]
+query query_boundary = query' []
+    where query' out q
+            | not $ query_boundary `intersects` (boundary_square q) = out
+            | otherwise = 
+                let es = filter (\e -> (boundary_square e) `intersects` query_boundary) $ elements q
+                in foldl (\out' cq -> query' out' cq) (out ++ es) (non_empty_children q)
+
diff --git a/src/System/Random/Utils.hs b/src/System/Random/Utils.hs
new file mode 100644
--- /dev/null
+++ b/src/System/Random/Utils.hs
@@ -0,0 +1,12 @@
+{-# LANGUAGE ScopedTypeVariables #-}
+module System.Random.Utils where
+import System.Random.Mersenne
+
+randomElement a gen = do
+    v :: Double <- random gen
+    let i :: Int = floor $ v * (fromIntegral $ length a)
+    return $ a !! i
+
+randomRange low high gen = do
+    v <- random gen
+    return $ low * (1 - v) + high * v
diff --git a/test/QuadTreeVisualize.hs b/test/QuadTreeVisualize.hs
new file mode 100644
--- /dev/null
+++ b/test/QuadTreeVisualize.hs
@@ -0,0 +1,101 @@
+{-# LANGUAGE ScopedTypeVariables #-}
+module Main where
+
+import Data.QuadTree
+
+import Data.AABB
+
+import Render ( init_display
+              , new_viewer
+              , view_rendering
+              )
+
+import Control.Monad
+
+import Data.Maybe ( maybe )
+
+import Graphics.Rendering.OpenGL as GL
+import Graphics.UI.GLUT
+
+import System.Random.Mersenne
+import System.Random.Utils
+
+main = do 
+    (viewer, _) <- new_viewer
+    gen <- newMTGen Nothing
+    (rq, gen) <- random_quadtree gen (empty :: QuadTree Elem) 5
+    view_quadtree viewer gen rq
+    return ()
+
+data Elem = Elem Boundary (Color3 Double)
+
+instance HasBoundary Elem where
+    boundary_points (Elem b _) = boundary_points b
+    boundary_edges (Elem b _) = boundary_edges b
+    boundary_extents (Elem b _) = boundary_extents b
+    boundary_square (Elem b _) = boundary_square b
+
+random_quadtree gen q 0 = return (q, gen)
+random_quadtree gen q n = do
+    x :: Double <- randomRange (-10.0) 10.0 gen
+    y :: Double <- randomRange (-10.0) 10.0 gen
+    s :: Double <- randomRange 0.001 1.0 gen
+    let eb = Boundary (x, y) s
+    r :: Double <- randomRange 0.0 1.0 gen
+    g :: Double <- randomRange 0.0 1.0 gen
+    b :: Double <- randomRange 0.0 1.0 gen
+    let e = Elem eb (Color3 r g b)
+    let q' = insert e q
+    random_quadtree gen q' (n - 1)
+    
+view_quadtree viewer gen q = do
+    view_rendering viewer $ display_quadtree gen q
+    mainLoop
+
+display_quadtree gen q@(QuadTree _ (Boundary (bx, by) bsize) _) = do
+    clearColor $= Color4 1.0 1.0 1.0 0.0
+    clear [ColorBuffer]
+    matrixMode $= Modelview 0
+    loadIdentity
+    -- First set the display to have the extents (0,0) and (1,1)
+    translate $ Vector3 (-1.0 :: Double) (-1.0) 0.0
+    scale (2.0 :: Double) 2.0 0.0
+    -- Scale the entire quadtree to the display.
+    scale (1.0 / bsize) (1.0 / bsize) 1.0
+    translate $ Vector3 (-bx) (-by) 0.0
+    lineWidth $= 2.0
+    render_elements q
+    outline_quadtree q
+    flush
+
+outline_quadtree q = do
+    color $ Color3 (0.0 :: Float) 0.0 0.0
+    polygonMode $= (Line, Line)
+    renderPrimitive Quads $ outline_quadtree' q
+    where 
+        outline_quadtree' (QuadTree _ b (cq0, cq1, cq2, cq3)) = do
+            mapM_ (maybe (return ()) $ \cq -> outline_quadtree' cq) [cq0, cq1, cq2, cq3]
+            render_boundary b
+
+render_elements q = do
+    polygonMode $= (Fill, Fill)
+    renderPrimitive Quads $ outline_elements' q
+    where
+        outline_elements' (QuadTree es _ (cq0, cq1, cq2, cq3)) = do
+            mapM_ (maybe (return ()) $ \cq -> outline_elements' cq) [cq0, cq1, cq2, cq3]
+            forM_ es $ \(Elem b c) -> do
+                color c
+                render_boundary b
+
+random_color_gen gen = sequence $ repeat $ do
+    r :: Double <- randomRange 0.0 1.0 gen
+    g :: Double <- randomRange 0.0 1.0 gen
+    b :: Double <- randomRange 0.0 1.0 gen
+    return $ Color3 r g b
+
+render_boundary (Boundary (x,y) size) = do
+    vertex $ Vertex2 x y
+    vertex $ Vertex2 (x + size) y
+    vertex $ Vertex2 (x + size) (y + size)
+    vertex $ Vertex2 x (y + size)
+
diff --git a/test/Render.hs b/test/Render.hs
new file mode 100644
--- /dev/null
+++ b/test/Render.hs
@@ -0,0 +1,163 @@
+{-# LANGUAGE NoMonomorphismRestriction, ScopedTypeVariables #-}
+module Render where
+
+import Control.Monad ( forM, when )
+
+import Data.Either
+import Data.Maybe
+
+import Foreign.Ptr
+
+import Data.IORef
+
+import Graphics.Rendering.OpenGL as GL
+import Graphics.UI.GLUT
+
+import System.Environment
+
+data Viewer = Viewer
+    {
+        image_dirty :: Bool,
+        mwindow_size :: Maybe Size,
+        mrendering :: Maybe (IO ()),
+        fbo_state :: FBOState,
+        viewing_position :: (Int, Int),
+        viewing_scale :: Double
+    }
+
+data FBOState = 
+      NoFBO
+    | HaveFBO FBOContext
+    deriving (Eq, Show)
+
+data FBOContext = FBOContext
+    {
+        fbo :: FramebufferObject,
+        fbo_tex :: TextureObject,
+        fbo_size :: Size
+    }
+    deriving (Eq, Show)
+
+init_display = do
+    name <- getProgName
+    args <- getArgs
+    args' <- initialize name args
+    initialDisplayMode $= [RGBAMode, WithDepthBuffer, DoubleBuffered]
+    createWindow "Render"
+    return args'
+
+new_viewer = do
+    args <- init_display
+    viewer_ref <- newIORef $ Viewer True Nothing Nothing NoFBO (0,0) 1.0
+    displayCallback $= viewer_display viewer_ref
+    reshapeCallback $= Just (reshape_viewer_window viewer_ref)
+    keyboardMouseCallback $= Just (glut_event_handler viewer_ref)
+    return (viewer_ref, args)
+
+reshape_viewer_window :: IORef Viewer -> ReshapeCallback
+reshape_viewer_window viewer_ref new_window_size = modifyIORef viewer_ref $ \viewer -> viewer { mwindow_size = Just new_window_size }
+
+view_rendering viewer_ref r = modifyIORef viewer_ref $ \viewer -> viewer { image_dirty = True, mrendering = Just r }
+
+initialize_fbo viewer_ref = do
+    [fbo :: FramebufferObject] <- genObjectNames 1
+    [fbo_tex :: TextureObject] <- genObjectNames 1
+    bindFramebufferEXT FramebufferTarget fbo
+    textureBinding Texture2D $= Just fbo_tex
+    texImage2D Nothing NoProxy 0 RGBA8 (TextureSize2D 4096 4096) 0 (PixelData RGBA UnsignedByte nullPtr)
+    textureFilter Texture2D $= ((Nearest, Nothing), Nearest)
+    framebufferTexture2DEXT FramebufferTarget (ColorAttachment 0) Texture2D fbo_tex 0
+    textureBinding Texture2D $= Nothing
+    let fbo_context = FBOContext fbo fbo_tex (Size 4096 4096)
+    modifyIORef viewer_ref $ \viewer -> viewer { fbo_state = HaveFBO fbo_context }
+    unbindFramebufferEXT FramebufferTarget
+
+rasterize_rendering viewer_ref = do
+    HaveFBO fbo_context <- readIORef viewer_ref >>= return . fbo_state
+    bindFramebufferEXT FramebufferTarget $ fbo fbo_context
+    preservingAttrib [ViewportAttributes] $ do
+        viewport $= (Position 0 0, fbo_size fbo_context)
+        fromJust . mrendering =<< readIORef viewer_ref
+        finish
+    unbindFramebufferEXT FramebufferTarget
+    modifyIORef viewer_ref $ \viewer -> viewer { image_dirty = False }
+    
+viewer_display viewer_ref = do
+    whenM (readIORef viewer_ref >>= return . isJust . mwindow_size) $ do
+        whenM (readIORef viewer_ref >>= return . (== NoFBO) . fbo_state) $ initialize_fbo viewer_ref
+        whenM ( do rerender <-  readIORef viewer_ref >>= return . image_dirty
+                   can_render <- readIORef viewer_ref >>= return . isJust . mrendering
+                   return $ can_render && rerender
+              ) $ rasterize_rendering viewer_ref
+        window_size <- return . fromJust . mwindow_size =<< readIORef viewer_ref
+        HaveFBO fbo_context <- readIORef viewer_ref >>= return . fbo_state
+        (pos_x, pos_y) <- readIORef viewer_ref >>= return . viewing_position
+        s <- readIORef viewer_ref >>= return . viewing_scale
+        viewport $= (Position 0 0, window_size)
+        -- Establish the coordinate center (0,0) -> (1.0, 1.0)
+        matrixMode $= Modelview 0
+        loadIdentity
+        scale s s 1
+        translate $ Vector3 (-1.0 :: Double) (-1.0) 0.0
+        scale (2.0 :: Double) 2.0 0.0
+        -- 
+        scale (1.0 / fromIntegral (width window_size)) 
+              (1.0 / fromIntegral (height window_size)) 
+              (1.0 :: Double)
+        translate $ Vector3 (-1.0 * fromIntegral pos_x) (-1.0 * fromIntegral pos_y) (0.0 :: Double)
+        -- Now render the FBO image.
+        clearColor $= Color4 1.0 1.0 1.0 0.0
+        clear [ColorBuffer]
+        texture Texture2D $= Enabled
+        textureBinding Texture2D $= Just (fbo_tex fbo_context)
+        textureFunction $= Replace
+        polygonMode $= (Fill, Fill)
+        let max_x :: Double = fromIntegral $ width $ fbo_size fbo_context
+        let max_y :: Double = fromIntegral $ height $ fbo_size fbo_context
+        renderPrimitive Quads $ do
+            color $ Color3 (0.0 :: Double) 0.0 0.0
+            vertex $ Vertex2 (0.0 :: Double) 0.0
+            texCoord $ TexCoord2 (0.0 :: Double) 0.0
+            vertex $ Vertex2 (0.0 :: Double) max_y
+            texCoord $ TexCoord2 (0.0 :: Double) 1.0
+            vertex $ Vertex2 max_x max_y
+            texCoord $ TexCoord2 (1.0 :: Double) 1.0
+            vertex $ Vertex2 max_x 0.0
+            texCoord $ TexCoord2 (1.0 :: Double) 0.0
+        textureBinding Texture2D $= Nothing
+        finish
+        swapBuffers
+        return ()
+
+glut_event_handler viewer_ref key key_state mods pos = do
+    (pos_x, pos_y) <- readIORef viewer_ref >>= return . viewing_position
+    let (pos_x', pos_y') = 
+            case (key,key_state,mods,pos) of
+                (SpecialKey KeyRight, Down, Modifiers Up Up Up, _) -> (pos_x + 10, pos_y)
+                (SpecialKey KeyLeft, Down, Modifiers Up Up Up, _) -> (pos_x - 10, pos_y)
+                (SpecialKey KeyUp, Down, Modifiers Up Up Up, _) -> (pos_x, pos_y + 10)
+                (SpecialKey KeyDown, Down, Modifiers Up Up Up, _) -> (pos_x, pos_y - 10)
+                otherwise -> (pos_x, pos_y)
+    s <- readIORef viewer_ref >>= return . viewing_scale
+    let s' = case (key,key_state,mods,pos) of
+                (SpecialKey KeyUp, Down, Modifiers Down Up Up, _) -> s * 2
+                (SpecialKey KeyDown, Down, Modifiers Down Up Up, _) -> s * 0.5
+                otherwise -> s
+    modifyIORef 
+        viewer_ref 
+        $ \viewer -> viewer 
+            { 
+                viewing_position = (pos_x', pos_y'),
+                viewing_scale = s'
+            }
+    postRedisplay Nothing
+
+whenM :: Monad m => m Bool -> m () -> m ()
+whenM bm m = bm >>= (\b -> when b m)
+
+width :: Size -> GLsizei
+width (Size w _) = w
+
+height :: Size -> GLsizei
+height (Size _ h) = h
+
diff --git a/test/Verify.hs b/test/Verify.hs
new file mode 100644
--- /dev/null
+++ b/test/Verify.hs
@@ -0,0 +1,14 @@
+module Main where
+
+import Verify.Data.QuadTree
+import Verify.Data.AABB
+import Test.QuickCheck
+
+main = do
+    quickCheck $ label "intersects_is_reflexive_prop" intersects_is_reflexive_prop
+    quickCheck $ label "encloses_is_reflexive_prop" encloses_is_reflexive_prop
+    quickCheck $ label "element_bounds_query_is_element_prop" element_bounds_query_is_element_prop
+    quickCheck $ label "oob_bounds_query_is_empty_prop" oob_bounds_query_is_empty_prop
+    quickCheck $ label "all_elements_inserted_query_prop " all_elements_inserted_query_prop 
+    putStrLn "DONE"
+
diff --git a/test/Verify/Data/AABB.hs b/test/Verify/Data/AABB.hs
new file mode 100644
--- /dev/null
+++ b/test/Verify/Data/AABB.hs
@@ -0,0 +1,45 @@
+module Verify.Data.AABB ( module Verify.Data.AABB
+                        , module Data.AABB
+                        )
+
+    where
+
+import Data.AABB
+
+import Control.Monad
+import Test.QuickCheck
+
+instance Arbitrary Boundary where
+    arbitrary = do
+        corner <- arbitrary
+        s <- liftM abs arbitrary
+        return $ Boundary corner s
+
+data NonIntersectingBounds = NonIntersectingBounds Boundary Boundary
+    deriving (Eq, Show)
+
+{- Two non-intersecting bounds can be generated by generating an arbitrary rectange defined by a min
+ - extent and max extent.
+ - The min extent is the corner of one boundary, p_0. The max extent is the corner of the other boundary, p_1.
+ - The boundary with a corner at p_1 can have any size.
+ - While the boundary with a corner at p_0 can not be given a size that could imply a boundary
+ - intersection.
+ -}
+instance Arbitrary NonIntersectingBounds where
+    arbitrary = do
+        x_0 <- arbitrary
+        y_0 <- arbitrary
+        x_1 <- arbitrary
+        y_1 <- arbitrary
+        let p_0 = (min x_0 x_1, min y_0 y_1)
+            p_1 = (max x_0 x_1, max y_0 y_1)
+        s_0 <- choose (0.0, min (fst p_1 - fst p_0) (snd p_1 - snd p_0))
+        s_1 <- arbitrary
+        return $ NonIntersectingBounds (Boundary p_0 s_0) (Boundary p_1 s_1)
+
+intersects_is_reflexive_prop :: Boundary -> Bool
+intersects_is_reflexive_prop b = b `intersects` b
+
+encloses_is_reflexive_prop :: Boundary -> Bool
+encloses_is_reflexive_prop b = b `encloses` b
+
diff --git a/test/Verify/Data/QuadTree.hs b/test/Verify/Data/QuadTree.hs
new file mode 100644
--- /dev/null
+++ b/test/Verify/Data/QuadTree.hs
@@ -0,0 +1,61 @@
+{-# LANGUAGE FlexibleInstances, ScopedTypeVariables #-}
+module Verify.Data.QuadTree ( module Data.QuadTree
+                            , module Verify.Data.QuadTree
+                            )
+    where
+
+import Data.QuadTree
+
+import Verify.Data.AABB
+
+import Control.Monad
+import Data.List hiding (insert)
+import Test.QuickCheck
+
+data BoundaryQTConstruction = BoundaryQTConstruction [Boundary] (QuadTree Boundary)
+    deriving (Show)
+
+instance Arbitrary BoundaryQTConstruction where
+    arbitrary = do
+        element_count <- choose (1,100)
+        (q, es) <- foldM (\(q, es) _ -> do
+                           e <- arbitrary
+                           return (insert e q, e : es)
+                         ) 
+                         (empty, []) 
+                         [1 :: Int .. element_count]
+        return $ BoundaryQTConstruction es q
+
+-- All elements inserted into the quadtree should be returned by a query for all elements within the
+-- boundaries of the quadtree
+all_elements_inserted_query_prop :: BoundaryQTConstruction -> Bool
+all_elements_inserted_query_prop (BoundaryQTConstruction es q) = es \\ query (boundary_square q) q == []
+
+element_bounds_query_is_element_prop :: Boundary -> Boundary -> Property
+element_bounds_query_is_element_prop initial_bounds element_bounds = 
+    boundary_size initial_bounds /= 0.0 ==>
+    let q = empty_with_bounds initial_bounds
+        q' = insert element_bounds q
+    in case query element_bounds q' of
+        []  -> False
+        [e] -> e == element_bounds
+        _   -> False
+
+oob_bounds_query_is_empty_prop :: NonIntersectingBounds -> Property
+oob_bounds_query_is_empty_prop (NonIntersectingBounds b_0 b_1) = 
+    boundary_size b_0 /= 0.0 ==>
+    let q :: QuadTree Boundary = empty_with_bounds b_0
+    in  [] ==  query b_1 q 
+
+-- An easy quadtree to test is one where the elements contained in the quadtree are boundaries.
+instance Show (QuadTree Boundary) where
+    show (QuadTree es b cq) = show es ++ " " ++ show b ++ " " ++ show cq ++ "\n"
+
+{-
+instance Show ( Maybe (QuadTree Boundary), Maybe (QuadTree Boundary)
+              , Maybe (QuadTree Boundary), Maybe (QuadTree Boundary) ) where
+    show (mq0, mq1, mq2, mq3) = "( " ++ show (fmap show mq0) ++ "," ++ show (fmap show mq1) ++
+                                "," ++ show (fmap show mq2) ++ "," ++ show (fmap show mq3) ++ ")"
+
+-}
+
diff --git a/test/run_verify b/test/run_verify
new file mode 100644
--- /dev/null
+++ b/test/run_verify
@@ -0,0 +1,7 @@
+#!/usr/bin/env runhaskell
+import System.Console.GetOpt
+import System.Cmd
+import System.Environment
+import Control.Monad
+
+main = system "ghc -ignore-package data-spacepart -i../src -o VerifyDriver --make Verify.hs && ./VerifyDriver"
diff --git a/test/run_visualize b/test/run_visualize
new file mode 100644
--- /dev/null
+++ b/test/run_visualize
@@ -0,0 +1,7 @@
+#!/usr/bin/env runhaskell
+import System.Console.GetOpt
+import System.Cmd
+import System.Environment
+import Control.Monad
+
+main = system "ghc -O3 --make QuadTreeVisualize.hs && ./QuadTreeVisualize"
