diff --git a/Data/Octree.hs b/Data/Octree.hs new file mode 100644 --- /dev/null +++ b/Data/Octree.hs @@ -0,0 +1,15 @@ +module Data.Octree(Octree, + Vector3(..), + dist, + fromList, toList, + lookup, + insert, + nearest, + depth, + size, + withinRange) +where + +import Prelude hiding(lookup) +import Data.Octree.Internal + diff --git a/Data/Octree/BoundingBox/BoundingBox.hs b/Data/Octree/BoundingBox/BoundingBox.hs new file mode 100644 --- /dev/null +++ b/Data/Octree/BoundingBox/BoundingBox.hs @@ -0,0 +1,67 @@ +{- | + Module : Data.Octree.BoundingBox.BoundingBox + Copyright : Copyright (c) 2016 Michael Litchard + License : BSD3 + + Maintainer : Michael Litchard + Stability : experimental + Portability: not portable + + This module provides a way to use bounding boxes with Octrees. + +-} + +module Data.Octree.BoundingBox.BoundingBox + ( BBoxConfig (..) + , traverseOctreeBB + , defBBoxConfig + ) where + +import Data.Maybe (mapMaybe) +import Data.BoundingBox.B3 (BBox3) + +import Data.Octree.Internal (Octree (..), allOctants) +import Data.Octree.BoundingBox.Internal + +-- | BBoxConfig - The functions traverseOctreeBB needs +data BBoxConfig x y a = BBoxConfig { +-- | A function to recurse down the Octree + select :: BBox3 -> x -> Maybe x, +-- | A function to pre-condition the leaves + procLeaf :: BBox3 -> x -> [LeafValue a] -> y, +-- | A function to recurse back up the tree + combine :: x -> [y] -> y +} + +-- | defBBoxConfig - default BBoxConfig +defBBoxConfig :: BBoxConfig DefInput DefOutput DefNodeValue +defBBoxConfig = BBoxConfig { + select = filterNodes , + procLeaf = points, + combine = result +} + +-- | traverseOctreeBB - Generalized Octree traversal function based on BBox3 +-- Base case processes leaves +traverseOctreeBB :: BBoxConfig x y a -> BBox3 -> Octree a -> x -> y +traverseOctreeBB bbc bbx (Leaf leaf_vals) input = procLeaf' bbx input leaf_vals + where + procLeaf' = procLeaf bbc + +-- | General framework to traverse an Octree in terms of a BBox3 +traverseOctreeBB + bbc bbx (Node split' nwu' nwd' neu' ned' swu' swd' seu' sed') x = + combine' x res + where + combine' = combine bbc + res = mapMaybe traverseOctreeBB' octList + select' = select bbc + + traverseOctreeBB' (subbox, subtree) = + case select' subbox x of + Just x' -> Just (traverseOctreeBB bbc subbox subtree x) + Nothing -> Nothing + + octList = zip boxList children + boxList = map (newBBox3 bbx split') allOctants + children = [swd',sed',nwd',ned',swu',seu',nwu',neu'] diff --git a/Data/Octree/BoundingBox/Internal.hs b/Data/Octree/BoundingBox/Internal.hs new file mode 100644 --- /dev/null +++ b/Data/Octree/BoundingBox/Internal.hs @@ -0,0 +1,122 @@ +{- | + Module : Data.Octree.BoundingBox.Internal + Copyright : Copyright (c) 2016 Michael Litchard + License : BSD3 + + Maintainer : Michael Litchard + Stability : experimental + Portability: not portable + + This module provides functions etc, for Data.Octree.BoundingBox.BoundingBox + + -- | DefInput, LeafValue a, DefOutput a, DefNodeValue a + Default types for BBoxConfig + + -- | filterNodes, points, result + Default functions for BBoxConfig + + -- | inclusive + boolean check for one BBox3 being inclusive of another +-} + +module Data.Octree.BoundingBox.Internal + ( filterNodes + , points + , result + , DefInput + , LeafValue + , DefOutput + , DefNodeValue + , newBBox3 + , inclusive + ) where + +import Data.BoundingBox.B3 (BBox3 (..), bound_corners, within_bounds, isect) + +import Data.Octree.Internal (Vector3 (..), ODir (..)) +import Data.List (foldl') + +type DefInput = Vector3 +type Split = Vector3 +type LeafValue a = (Vector3, a) +type DefOutput = (BBox3, [LeafValue DefNodeValue ]) +type DefNodeValue = Int + +-- | newBBox3 creates a smaller BBox3 +-- Given Node name, previous BBox3, and the split +newBBox3 :: BBox3 -> Split -> ODir -> BBox3 +newBBox3 bbx split' SWD = + bound_corners swdCorner neuCorner + where + swdCorner = Vector3 (minX bbx) (minY bbx) (minZ bbx) + neuCorner = Vector3 (v3x split') (v3y split') (v3z split') + +newBBox3 bbx split' SED = + bound_corners swdCorner neuCorner + where + swdCorner = Vector3 (v3x split') (minY bbx) (minZ bbx) + neuCorner = Vector3 (maxX bbx) (v3y split') (v3z split') + +newBBox3 bbx split' NWD = + bound_corners swdCorner neuCorner + where + swdCorner = Vector3 (minX bbx) (v3y split') (minZ bbx) + neuCorner = Vector3 (v3x split') (maxY bbx) (v3z split') + +newBBox3 bbx split' NED = + bound_corners swdCorner neuCorner + where + swdCorner = Vector3 (v3x split') (v3y split') (minZ bbx) + neuCorner = Vector3 (maxX bbx) (maxY bbx) (v3z split') + +newBBox3 bbx split' SWU = + bound_corners swdCorner neuCorner + where + swdCorner = Vector3 (minX bbx) (minY bbx) (v3z split') + neuCorner = Vector3 (v3x split') (v3y split') (maxZ bbx) + +newBBox3 bbx split' SEU = + bound_corners swdCorner neuCorner + where + swdCorner = Vector3 (v3x split') (minY bbx) (v3z split') + neuCorner = Vector3 (maxX bbx) (v3y split') (maxZ bbx) + +newBBox3 bbx split' NWU = + bound_corners swdCorner neuCorner + where + swdCorner = Vector3 (minX bbx) (v3y split') (v3z split') + neuCorner = Vector3 (v3x split') (maxY bbx) (maxZ bbx) + +newBBox3 bbx split' NEU = + bound_corners swdCorner neuCorner + where + swdCorner = Vector3 (v3x split') (v3y split') (v3z split') + neuCorner = Vector3 (maxX bbx) (maxY bbx) (maxZ bbx) + +-- | filterNodes is default function for BBoxConfig +-- used to recurse down octree identifying which BBox3s contain DefInput +filterNodes :: BBox3 -> DefInput -> Maybe DefInput +filterNodes bbox x = if within_bounds x bbox then Just x else Nothing + +-- | points is default function for BBoxConfig +-- pre-processes Leaf +points :: BBox3 -> DefInput -> [LeafValue DefNodeValue ] -> DefOutput +points box _ leaf = (box, leaf) + +-- | result reduces the list of all BBoxes containing the point to +-- the terminal BBox +result :: DefInput -> [DefOutput ] -> DefOutput +result _ (x:xs) = foldl' findTerminal x xs + where + findTerminal :: DefOutput -> DefOutput -> DefOutput + findTerminal bbox1@(bbox1',_) bbox2@(bbox2',_) + | inclusive bbox1' bbox2' = bbox1 + | otherwise = bbox2 + +-- | Supplied boolean test for result function +-- Returns True if Box b1 is contained within Box b2 +inclusive :: BBox3 -> BBox3 -> Bool +inclusive b1 b2 = + case isect b1 b2 of + Just b3 -> b1 == b3 + Nothing -> False diff --git a/Data/Octree/BoundingBox/README.lhs b/Data/Octree/BoundingBox/README.lhs new file mode 100644 --- /dev/null +++ b/Data/Octree/BoundingBox/README.lhs @@ -0,0 +1,53 @@ +point-octree/Data/Octree/BoundingBox +============================== +BoundingBox.hs provides a generalized traversal function, that traverses +an `Octree a` based on a BBox3 value. + +Here is an example that produces an `Octree Int` with 2 subdivisions, +and finds the bounding box of a particular `(Vector3,Int)`. + +~~~ {.haskell} +module Main where + +import Data.List +import Data.Maybe +import Data.Vector.Class +import System.Random +import System.Random.Shuffle +import Data.BoundingBox.B3 + +import Data.Octree.Internal hiding (lookup) +import Data.Octree.BoundingBox.BoundingBox +import Data.Octree.BoundingBox.Internal + +main = do + octree' <- octree 513 + let (testInput,_) = (!!) (toList octree') 256 + (bbx,_) = traverseOctreeBB defBBoxConfig bbis octree' testInput + putStrLn ("Bounding box of " ++ (show testInput) ++ " is " ++ (show bbx)) + +bbis :: BBox3 +bbis = + let infinity = (read "Infinity") :: Double + swdCorner = Vector3 (-infinity) (-infinity) (-infinity) + neuCorner = Vector3 (infinity) (infinity) (infinity) + in bound_corners swdCorner neuCorner + +octree :: Integer -> IO (Octree DefNodeValue) +octree bound = do + xGen <- newStdGen + yGen <- newStdGen + zGen <- newStdGen + let xPoints :: [Double] + yPoints :: [Double] + zPoints :: [Double] + xPoints = map fromInteger $ shuffle' pointList (length pointList) xGen + yPoints = map fromInteger $ shuffle' pointList (length pointList) yGen + zPoints = map fromInteger $ shuffle' pointList (length pointList) zGen + ub = bound `div` 2 + lb = - (bound `div` 2) + pointList = [lb .. ub] + sizePL = length pointList + vectors = map (\(x,y,z) -> Vector3 x y z) $ zip3 xPoints yPoints zPoints + return $ fromList $ zip vectors $ [1 .. sizePL] +~~~ diff --git a/Data/Octree/Internal.hs b/Data/Octree/Internal.hs new file mode 100644 --- /dev/null +++ b/Data/Octree/Internal.hs @@ -0,0 +1,356 @@ +{- | + Module : Data.Octree.Internal + Copyright : Copyright (c) 2012 Michal J. Gajda + License : BSD3 + + Maintainer : Michael Litchard + Stability : experimental + Portability : not portable + + This module provides the description of a point octree + and it's operations. + +-} +{-# LANGUAGE ScopedTypeVariables, DisambiguateRecordFields #-} +{-# LANGUAGE DeriveFunctor, DeriveFoldable, DeriveTraversable #-} +module Data.Octree.Internal + ( Vector3(..), dist + , Octree(..), lookup, nearest, withinRange, fromList, toList, insert + -- internal + , ODir(..) + , octreeStep, octantDistance, splitBy', joinStep, splitStep, allOctants, octantDistance' + , cmp + , pickClosest + , depth, size + , subnodes + ) where + +import Data.Functor (Functor (..)) +import Data.Foldable (Foldable (foldr)) +import Data.Traversable(Traversable(..)) + +import Data.Vector.V3 +import Data.Vector.Class + +import Prelude hiding(lookup, foldr) +import Data.List(sort, sortBy) +import Data.Maybe(maybeToList, listToMaybe) +import Data.Bits((.&.)) +import Control.Arrow(second) +import Test.QuickCheck.All(quickCheckAll) +import Test.QuickCheck.Arbitrary + +-- | norm of a vector +norm :: Vector3 -> Double +norm a = sqrt (a `vdot` a) + +-- | distance between two vectors +dist :: Vector3 -> Vector3 -> Double +dist u v = norm (u - v) + +-- | Datatype for nodes within Octree. +data Octree a + = Node { split :: Vector3 + , nwu :: Octree a + , nwd :: Octree a + , neu :: Octree a + , ned :: Octree a + , swu :: Octree a + , swd :: Octree a + , seu :: Octree a + , sed :: Octree a + } + | Leaf { unLeaf :: [(Vector3, a)] } + deriving (Show, Functor, Foldable, Traversable) + +-- | Enumerated type to indicate octants in 3D-space relative to given center. +data ODir + = SWD + | SED + | NWD + | NED + | SWU + | SEU + | NWU + | NEU + deriving (Eq, Ord, Enum, Show, Bounded) + +-- | Convinience synonym to make function signatures more readable +type OctPoints a = ( [(Vector3, a)] + , [(Vector3, a)] + , [(Vector3, a)] + , [(Vector3, a)] + , [(Vector3, a)] + , [(Vector3, a)] + , [(Vector3, a)] + , [(Vector3, a)]) + +-- | Internal method that gives octant of a first vector, +-- relative to the second vector as a center. +cmp :: Vector3 -> Vector3 -> ODir +cmp ca cb = joinStep (cx, cy, cz) + where + cx = v3x ca >= v3x cb + cy = v3y ca >= v3y cb + cz = v3z ca >= v3z cb + +-- | Internal method that joins result of three coordinate comparisons, +-- and makes an octant name `ODir` +joinStep :: (Enum a1, Enum a3, Enum a2, Enum a) => (a1, a2, a3) -> a +joinStep (cx, cy, cz) = toEnum (fromEnum cx + 2 * fromEnum cy + 4 * fromEnum cz) + +-- | This function converts octant name to a function +-- that steps down in an Octree towards this octant +octreeStep :: Octree a -> ODir -> Octree a +octreeStep ot NWU = nwu ot +octreeStep ot NWD = nwd ot +octreeStep ot NEU = neu ot +octreeStep ot NED = ned ot +octreeStep ot SWU = swu ot +octreeStep ot SWD = swd ot +octreeStep ot SEU = seu ot +octreeStep ot SED = sed ot + +-- | Function that splits octant name into three boolean values, +-- depending of sign of a relative coordinate in that octant. +-- (Coordinate is relative to a split point within Octree.) +splitStep :: ODir -> (Bool, Bool, Bool) +splitStep step = ((val .&. 1) == 1, (val .&. 2) == 2, (val .&. 4) == 4) + where val = fromEnum step + +-- | Internal function that finds a lower bound for a distance +-- between a point of relative coordinates, and octant of given name. +-- It works only when coordinates of a given point are all positive. +-- It is used only by `octantDistance`, which respectively changes +-- octant name depending of signs of relative coordinates. +-- Here we assume that a, b, c > 0 +-- (otherwise we will take abs, and correspondingly invert results) +-- same octant +-- dp = difference between given point and the center of Octree node +octantDistance' :: Vector3 -> ODir -> Scalar +octantDistance' dp NEU = 0.0 +-- | adjacent by plane +octantDistance' dp NWU = v3x dp +octantDistance' dp SEU = v3y dp +octantDistance' dp NED = v3z dp +-- | adjacent by edge +octantDistance' dp SWU = sqrt ( v3x dp * v3x dp + v3y dp * v3y dp) +octantDistance' dp SED = sqrt ( v3y dp * v3y dp + v3z dp * v3z dp) +octantDistance' dp NWD = sqrt ( v3x dp * v3x dp + v3z dp * v3z dp) +-- adjacent by point +octantDistance' dp SWD = norm dp + +-- | List of all octant names. +allOctants :: [ODir] +allOctants = [minBound..maxBound] + +-- | Internal function that makes code clearer. +xor :: Bool -> Bool -> Bool +xor = (/=) + +-- | Finds a minimum bounds for a distance between a given point +-- | in relative coordinates and a given octant. +octantDistance :: Vector3 -> ODir -> Double +octantDistance dp odir = octantDistance' (abs dp) (toggle dp odir) + +-- | Toggles octant names depending on a signs of vector coordinates +-- | for use in octantDistance. +toggle :: Vector3 -> ODir -> ODir +toggle dp odir = + joinStep ((v3x dp >= 0) `xor` not u, + (v3y dp >= 0) `xor` not v, + (v3z dp >= 0) `xor` not w) + where (u, v, w) = splitStep odir + +-- | Given a point in relative coordinates, gives list of all octants +-- and minimum distances from this point. +octantDistances :: Vector3 -> [(ODir, Double)] +octantDistances dp = [(o, octantDistance dp o) | o <- allOctants] + +-- | splits a list of vectors and "payload" tuples +-- | into a tuple with elements destined for different octants. +-- FIXME: VERY IMPORTANT - add prop_splitBy vs cmp +splitBy :: Vector3 -> [(Vector3, a)] -> OctPoints a +splitBy _splitPoint [] = ([], [], [], [], [], [], [], []) +splitBy splitPoint ((pt@(coord, a)):aList) = case i of + SWD -> (pt:swd, sed, nwd, ned, swu, seu, nwu, neu) + SED -> ( swd, pt:sed, nwd, ned, swu, seu, nwu, neu) + NWD -> ( swd, sed, pt:nwd, ned, swu, seu, nwu, neu) + NED -> ( swd, sed, nwd, pt:ned, swu, seu, nwu, neu) + SWU -> ( swd, sed, nwd, ned, pt:swu, seu, nwu, neu) + SEU -> ( swd, sed, nwd, ned, swu, pt:seu, nwu, neu) + NWU -> ( swd, sed, nwd, ned, swu, seu, pt:nwu, neu) + NEU -> ( swd, sed, nwd, ned, swu, seu, nwu, pt:neu) + where + i = cmp coord splitPoint + (swd, sed, nwd, ned, swu, seu, nwu, neu) = splitBy splitPoint aList + +-- | Computes a center of mass for a given list of vectors. +-- Used to find a splitPoint. +massCenter :: Fractional a => [(a, b)] -> a +massCenter aList = sum (map fst aList) / count + where count = fromInteger . toInteger . length $ aList + +-- | Helper function to map over an 8-element tuple +tmap :: (t -> t1) -> + (t, t, t, t, t, t, t, t) -> + (t1, t1, t1, t1, t1, t1, t1, t1) +tmap t (a, b, c, d, e, f, g, h) = (t a, t b, t c, t d, t e, t f, t g, t h) + +-- | Maximum number of elements before Octree leaf is split. +leafLimit :: Int +leafLimit = 16 + +-- | Creates an Octree from a list of (index, payload) tuples. +fromList :: [(Vector3, a)] -> Octree a +fromList aList = + if length aList <= leafLimit + then Leaf aList + else let splitPoint :: Vector3 = massCenter aList + in splitBy' fromList splitPoint aList + +-- | Internal method, that splits a list into octants depending on coordinates, +-- and then applies a specified function to each of these sublists, +-- in order to create subnodes of the Octree +splitBy' :: ([(Vector3, a)] -> Octree a1)-> Vector3-> [(Vector3, a)]-> Octree a1 +splitBy' f splitPoint aList = + Node { split = splitPoint + , nwu = tnwu + , nwd = tnwd + , neu = tneu + , ned = tned + , swu = tswu + , swd = tswd + , seu = tseu + , sed = tsed } + where + (tswd, tsed, tnwd, tned, tswu, tseu, tnwu, tneu) = + tmap f $ splitBy splitPoint aList +-- TODO: use arrays for memory savings + +-- | Internal method that prepends contents of the given subtree to a list +-- | given as argument. +toList' :: Octree t -> [(Vector3, t)] -> [(Vector3, t)] +toList' (Leaf l ) tmp = l ++ tmp +toList' (Node { nwu = a, + nwd = b, + neu = c, + ned = d, + swu = e, + swd = f, + seu = g, + sed = h }) tmp = foldr toList' tmp [a, b, c, d, e, f, g, h] +-- | Creates an Octree from list, trying to keep split points near centers +-- of mass for each subtree. +toList :: Octree t -> [(Vector3, t)] +toList t = toList' t [] + +-- | Finds a path to a Leaf where a given point should be, +-- and returns a list of octant names. +pathTo :: Vector3 -> Octree a -> [ODir] +pathTo pt (Leaf _) = [] +pathTo pt node = aStep : pathTo pt (octreeStep node aStep) + where aStep = cmp pt (split node) + +-- | Applies a given function to a node specified by a path +-- (list of octant names), and then returns a modified Octree. +applyByPath :: (Octree a -> Octree a) -> [ODir] -> Octree a -> Octree a +applyByPath f [] ot = f ot +applyByPath f (step:path) node = case step of + NWU -> node{ nwu = applyByPath f path (nwu node) } + NWD -> node{ nwd = applyByPath f path (nwd node) } + NEU -> node{ neu = applyByPath f path (neu node) } + NED -> node{ ned = applyByPath f path (ned node) } + SWU -> node{ swu = applyByPath f path (swu node) } + SWD -> node{ swd = applyByPath f path (swd node) } + SEU -> node{ seu = applyByPath f path (seu node) } + SED -> node{ sed = applyByPath f path (sed node) } + +-- | Inserts a point into an Octree. +-- NOTE: insert accepts duplicate points, +-- but lookup would not find them. Use withinRange in such case. +insert :: (Vector3, a) -> Octree a -> Octree a +insert (pt, dat) ot = applyByPath insert' path ot + where + path = pathTo pt ot + insert' (Leaf l) = fromList ((pt, dat) : l) + insert' _ = error "Impossible in insert'" + + +-- | Internal: finds candidates for nearest neighbour lazily for each octant; +-- they are returned in a list of (octant, min. bound for distance, +-- Maybe candidate) tuples. +candidates' :: Vector3 -> Octree a -> [(ODir, Double, [(Vector3, a)])] +candidates' pt (Leaf l) = [] +candidates' pt node = + map findCandidates . sortBy compareDistance . octantDistances $ + pt - split node + where + findCandidates (octant, d) = + (octant, d, maybeToList . pickClosest pt . toList . octreeStep node $ octant) + compareDistance a b = compare (snd a) (snd b) + +-- | Finds a given point, if it is in the tree. +lookup :: Octree a -> Vector3 -> Maybe (Vector3, a) +lookup (Leaf l) pt = listToMaybe . filter ((==pt) . fst) $ l +lookup node pt = flip lookup pt . octreeStep node . cmp pt . split $ node + +-- | Finds nearest neighbour for a given point. +nearest :: Octree a -> Vector3 -> Maybe (Vector3, a) +nearest (Leaf l) pt = pickClosest pt l +nearest node pt = selectFrom candidates + where + candidates = + map findCandidate . sortBy compareDistance . octantDistances $ + pt - split node + compareDistance a b = compare (snd a) (snd b) + findCandidate (octant, d) = + (maybeToList . nearest' . octreeStep node $ octant, d) + + nearest' n = nearest n pt + + selectFrom (([], _d) : cs) = selectFrom cs + selectFrom (([best], _d) : cs) = selectFrom' best cs + selectFrom [] = Nothing + selectFrom' best (([], d) : cs) = selectFrom' best cs + -- TODO: FAILS: shortcut guard to avoid recursion over + -- whole structure (since d is bound for distance within octant): + selectFrom' best ((c, d) : cs) | d > dist pt (fst best) = Just best + selectFrom' best (([next], d) : cs) = + selectFrom' nextBest cs + where nextBest = if dist pt (fst best) <= dist pt (fst next) + then best + else next + selectFrom' best [] = Just best + +-- | Internal method that picks from a given list a point closest to argument, +pickClosest :: Vector3 -> [(Vector3, t)] -> Maybe (Vector3, t) +pickClosest pt [] = Nothing +pickClosest pt (a:as) = Just $ foldr (pickCloser pt) a as +pickCloser pt va@(a, _a) vb@(b, _b) = if dist pt a <= dist pt b + then va + else vb + +-- | Returns all points within Octree that +-- are within a given distance from argument. +withinRange :: Octree a -> Scalar -> Vector3 -> [(Vector3, a)] +withinRange (Leaf l) r pt = filter (\(lpt, _) -> dist pt lpt <= r) l +withinRange node r pt = + concatMap recurseOctant . -- recurse over remaining octants, + -- and merge results + filter ((<=r) . snd) . -- discard octants that are out of range + octantDistances $ pt - split node -- find octant distances + where + recurseOctant (octant, _d) = + (\o -> withinRange o r pt) . octreeStep node $ octant + +subnodes :: Octree a -> [Octree a] +subnodes (Leaf _) = [] +subnodes node = map (octreeStep node) allOctants + +depth :: Octree a -> Int +depth (Leaf _) = 0 +depth node = foldr max 0 . map ((+1) . depth) . subnodes $ node + +size :: Octree a -> Int +size = length . toList diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,31 @@ +Haskell code in this package is subject to: + +Copyright (c) Michal J. Gajda '2012 + +Copyright (c) Michael A. Litchard '2016 + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions +are met: +1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. +2. 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. +3. Neither the name of the author nor the names of his contributors + may be used to endorse or promote products derived from this software + without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE AUTHOR 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 AUTHORS 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/README.lhs b/README.lhs new file mode 100644 --- /dev/null +++ b/README.lhs @@ -0,0 +1,26 @@ +point-octree +====== +This is a fork of Michał J Gajda's [Octree](https://github.com/BioHaskell/octree) library, with bounding boxes. + +[![Build Status](https://travis-ci.org/mlitchard/point-octree.svg?branch=master)](https://www.travis-ci.org/mlitchard/point-octree) + +(From Michal's README.md) +To use simply: + +~~~ {.haskell} +module Main where + +import Data.Octree as O + +import Data.Vector.V3 + +main = do let oct = fromList [(Vector3 1 2 3, "a"), + (Vector3 3 4 5, "b"), + (Vector3 8 8 8, "c")] + report msg elt = putStrLn $ msg ++ show elt + report "Nearest :" $ O.nearest oct $ Vector3 2 2 3 + report "Within range:" $ O.withinRange oct 5.0 $ Vector3 2 2 3 + return () +~~~ + + diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,4 @@ +#! /usr/bin/env runhaskell + +import Distribution.Simple +main = defaultMain diff --git a/changelog b/changelog new file mode 100644 --- /dev/null +++ b/changelog @@ -0,0 +1,3 @@ +-*-Changelog-*- +0.5.5.3 2016/07/03 + Hard fork from octree library diff --git a/point-octree.cabal b/point-octree.cabal new file mode 100644 --- /dev/null +++ b/point-octree.cabal @@ -0,0 +1,85 @@ +name: point-octree +version: 0.5.5.3 +stability: beta +homepage: https://github.com/mlitchard/point-octree +package-url: http://hackage.haskell.org/package/point-octree +synopsis: Point octree, with bounding boxes +description: Based on Michal J. Gajda's octree package, but with bounding boxes. +category: Data +license: BSD3 +license-file: LICENSE +extra-source-files: changelog + +author: Michal J. Gajda, Michael Litchard +copyright: Copyright by Michal J. Gajda '2012, Copyright by Michael Litchard '2016 +maintainer: michael@schmong.org +bug-reports: mailto:michael@schmong.org + + +build-type: Simple +cabal-version: >=1.8 +tested-with: GHC==7.6.3,GHC==7.8.3 + +source-repository head + type: git + location: git@github.com:mgajda/octree.git + +Library + build-depends: base >=4.0 && < 4.10, + AC-Vector >= 2.3.0, + QuickCheck >= 2.4.0 + + exposed-modules: Data.Octree + other-modules: Data.Octree.Internal + Data.Octree.BoundingBox.BoundingBox + Data.Octree.BoundingBox.Internal + exposed: True + extensions: ScopedTypeVariables + +Test-suite test_Octree + Type: exitcode-stdio-1.0 + build-depends: base >=4.0 && < 4.10, + hspec == 2.2.3, + AC-Vector >= 2.3.0, + QuickCheck >= 2.4.0 + + hs-source-dirs: ., Data, tests + Main-is: Main.hs + other-modules: PropTests.Octree + Data.Octree + Data.Octree.Internal + Data.Octree.BoundingBox.BoundingBox + Data.Octree.BoundingBox.Internal + PropTests.Common + PropTests.OctreeTests.Exposed + PropTests.OctreeTests.Internal + PropTests.BoundingBoxTests.BoundingBoxTests + PropTests.BoundingBoxTests.Utilities +Test-suite readme + type: exitcode-stdio-1.0 + -- We have a symlink: README.lhs -> README.md + main-is: README.lhs + Build-depends: base >=4.0 && < 4.10, + AC-Vector >= 2.3.0, + QuickCheck >= 2.4.0, + markdown-unlit + ghc-options: -pgmL markdown-unlit + other-modules: Data.Octree + Data.Octree.Internal + +Test-suite readme-BB + type: exitcode-stdio-1.0 + hs-source-dirs: . + main-is: Data/Octree/BoundingBox/README.lhs + Build-depends: base >=4.0 && < 4.10, + AC-Vector >= 2.3.0, + QuickCheck >= 2.4.0, + random == 1.1, + random-shuffle, + markdown-unlit + + ghc-options: -pgmL markdown-unlit + other-modules: Data.Octree + Data.Octree.Internal + Data.Octree.BoundingBox.BoundingBox + Data.Octree.BoundingBox.Internal diff --git a/tests/Main.hs b/tests/Main.hs new file mode 100644 --- /dev/null +++ b/tests/Main.hs @@ -0,0 +1,21 @@ +{- | + Module: tests/Main.hs + Copyright: 2016 Michael Litchard + License: BSD3 + Maintainer: + + Top-level test driver + +-} + +import Test.Hspec (hspec) + +import PropTests.Octree +import PropTests.BoundingBoxTests.BoundingBoxTests + +main :: IO () +main = do + hspec propOctExposed + hspec propValidateBoxes + hspec propMatchValuesToBoxes + hspec propOctInternal diff --git a/tests/PropTests/BoundingBoxTests/BoundingBoxTests.hs b/tests/PropTests/BoundingBoxTests/BoundingBoxTests.hs new file mode 100644 --- /dev/null +++ b/tests/PropTests/BoundingBoxTests/BoundingBoxTests.hs @@ -0,0 +1,72 @@ +{- | + Module : PropTests.BoundingBoxTests.BoundingBoxTests + Copyright : Copyright (c) 2016 Michael Litchard + License : BSD3 + + Maintainer : Michael Litchard + Stability : experimental + Portability: not portable + + This module provides property tests for Data.Octree.BoundingBox.BoundingBox + +-} + +module PropTests.BoundingBoxTests.BoundingBoxTests + ( propValidateBoxes + , propMatchValuesToBoxes + ) where + + + +import Test.Hspec (Spec, describe) +import Test.Hspec.QuickCheck (prop) + +import Data.Octree.Internal (Vector3 (..), fromList) +import Data.Octree.BoundingBox.BoundingBox (BBoxConfig, traverseOctreeBB) + +import PropTests.Common +import PropTests.BoundingBoxTests.Utilities +import Data.Octree.BoundingBox.Internal (DefInput, DefOutput, DefNodeValue) + +propValidateBoxes :: Spec +propValidateBoxes = + describe "prop_validateBoxes" $ + prop propMsg prop_validateBoxes + where + propMsg = + "Testing that arbitarily generated bounding boxes are well-formed." + +propMatchValuesToBoxes :: Spec +propMatchValuesToBoxes = + describe "prop_matchValuesToBox" $ + prop propMsg $ + prop_matchValuesToBox testConfig1 testConfig2 + where + propMsg = + "Testing that each generated BBox3 maps to correct [(Vector3, Int)]" + + +-- | validateBoxes takes the list of all bounding boxes of an Octree +-- and determines if they are all well-formed +prop_validateBoxes :: [(Vector3, Int)] -> Bool +prop_validateBoxes l = + let octree' = fromList l + bbxs = map fst $ traverseOctreeBB testConfig1 bbis octree' (Vector3 0 0 0) + in all isBBox bbxs + +type TestConfigOne = BBoxConfig DefInput [DefOutput] DefNodeValue +type TestConfigTwo = BBoxConfig AltInput DefOutput DefNodeValue +-- | matchValuesToBox verifies that each BBox3 generated by +-- TestConfigOne maps to the correct key list +prop_matchValuesToBox :: TestConfigOne -> + TestConfigTwo -> + [(Vector3, Int)] -> + Bool +prop_matchValuesToBox config1 config2 l = + let octree' = fromList l + boxCol1 = traverseOctreeBB config1 bbis octree' x + x = Vector3 0 0 0 -- this is never used but needs to be present + xs' = map fst boxCol1 + boxCol2 = map (traverseOctreeBB config2 bbis octree') xs' + validate = map (`elem` boxCol1) boxCol2 + in and validate diff --git a/tests/PropTests/BoundingBoxTests/Utilities.hs b/tests/PropTests/BoundingBoxTests/Utilities.hs new file mode 100644 --- /dev/null +++ b/tests/PropTests/BoundingBoxTests/Utilities.hs @@ -0,0 +1,88 @@ +{- | + Module : PropTests.BoundingBoxTests.Utilities + Copyright : Copyright (c) 2016 Michael Litchard + License : BSD3 + + Maintainer : Michael Litchard + Stability : experimental + Portability: not portable + + This module provides utility functions for + PropTests.BoundingBoxTests.External.hs + +-} + +module PropTests.BoundingBoxTests.Utilities + ( bbis + , isBBox + , testConfig1 + , testConfig2 + , AltInput + ) where + + +import Data.BoundingBox.B3 (BBox3 (..),bound_corners) + +import Data.Octree.Internal hiding (lookup) +import Data.Octree.BoundingBox.BoundingBox (BBoxConfig (..)) +import Data.Octree.BoundingBox.Internal ( DefInput + , DefOutput + , LeafValue + , DefNodeValue + , inclusive) +-- | Bounding Box of Infinite Space +-- The root Bounding Box +bbis :: BBox3 +bbis = + bound_corners swdCorner neuCorner + where + infinity = read "Infinity" :: Double + swdCorner = Vector3 (-infinity) (-infinity) (-infinity) + neuCorner = Vector3 infinity infinity infinity + +isBBox :: BBox3 -> Bool +isBBox (BBox3 minX minY minZ maxX maxY maxZ) = + (minX < maxX) && (minY < maxY) && (minZ < maxZ) + + +testConfig1 :: BBoxConfig DefInput [DefOutput] DefNodeValue +testConfig1 = BBoxConfig { + select = allPoints, + procLeaf = allBoxesAndKeys, + combine = assemble +} + +allBoxesAndKeys :: BBox3 -> + DefInput -> + [LeafValue DefNodeValue ] -> + [DefOutput] +allBoxesAndKeys bbx x leaf = [(bbx, leaf)] + +allPoints :: BBox3 -> DefInput -> Maybe DefInput +allPoints bbx = Just + +assemble :: DefInput -> [[DefOutput]] -> [DefOutput] +assemble _ = concat + +type AltInput = BBox3 + +testConfig2 :: BBoxConfig AltInput DefOutput DefNodeValue +testConfig2 = BBoxConfig { + select = filterBoxes, + procLeaf = checkbox, + combine = boxResult +} + +filterBoxes :: BBox3 -> AltInput -> Maybe AltInput +filterBoxes genBox origBox = + if inclusive origBox genBox then Just origBox else Nothing + +checkbox :: BBox3 -> AltInput -> [LeafValue DefNodeValue ] -> DefOutput +checkbox genBox _ leaf = (genBox, leaf) + +boxResult :: AltInput -> [DefOutput] -> DefOutput +boxResult origBox allOuts = + maybe bogusDummy ((,) origBox) $ lookup origBox allOuts + where + bogusDummy :: DefOutput + bogusDummy = (bbis,[]) diff --git a/tests/PropTests/Common.hs b/tests/PropTests/Common.hs new file mode 100644 --- /dev/null +++ b/tests/PropTests/Common.hs @@ -0,0 +1,34 @@ +{- | + Module : PropTests.Common + Copyright : Copyright (c) 2012 Michal J. Gajda + License : BSD3 + + Maintainer : Michael Litchard + Stability : experimental + Portability : not portable + + This module provides some helpers for the tests. + +-} + +module PropTests.Common where + +import Data.Vector.Class (vunpack) +import Data.Vector.V3 (Vector3 (Vector3)) +import Test.QuickCheck.Arbitrary + +-- | For testing purposes +instance Ord Vector3 where + a `compare` b = pointwiseOrd $ zipWith compare (vunpack a) (vunpack b) + +pointwiseOrd [] = EQ +pointwiseOrd (LT:cs) = LT +pointwiseOrd (GT:cs) = GT +pointwiseOrd (EQ:cs)= pointwiseOrd cs + +instance Arbitrary Vector3 where + arbitrary = do x <- arbitrary + y <- arbitrary + z <- arbitrary + return $ Vector3 x y z + diff --git a/tests/PropTests/Octree.hs b/tests/PropTests/Octree.hs new file mode 100644 --- /dev/null +++ b/tests/PropTests/Octree.hs @@ -0,0 +1,62 @@ +{- | + Module : PropTests.Octree + Copyright : Copyright (c) 2016 Michael Litchard + License : BSD3 + + Maintainer : Michael Litchard + Stability : experimental + Portability: not portable + + This module organizes the property tests for Octree proper. +-} + +{-# LANGUAGE ScopedTypeVariables #-} + +module PropTests.Octree + ( propOctInternal + , propOctExposed) + where + +import Data.Octree.Internal +import Data.Octree() -- test that interface module is not broken +import Prelude hiding(lookup) +import Data.List(sort, sortBy) + +import Test.Hspec (Spec, describe) +import Test.Hspec.QuickCheck (prop) + +import PropTests.OctreeTests.Internal +import PropTests.OctreeTests.Exposed + +propOctInternal :: Spec +propOctInternal = + describe "Tests for internal helper functions" $ do + prop "depth" prop_depth + prop "prop_cmp1" prop_cmp1 + prop "prop_cmp2" prop_cmp2 + prop "prop_stepDescription" prop_stepDescription + prop "prop_octantDistanceNoGreaterThanInterpointDistance0" + prop_octantDistanceNoGreaterThanInterpointDistance0 + prop "prop_octantDistanceNoGreaterThanInterpointDistance" + prop_octantDistanceNoGreaterThanInterpointDistance + prop "prop_octantDistanceNoGreaterThanInterpointDistanceZero" + prop_octantDistanceNoGreaterThanInterpointDistanceZero + prop "prop_octantDistanceNoGreaterThanCentroidDistance" + prop_octantDistanceNoGreaterThanCentroidDistance + prop "prop_pickClosest" + (prop_pickClosest :: [(Vector3, Int)] -> Vector3 -> Bool) + +propOctExposed :: Spec +propOctExposed = + describe "Tests for exposed functions" $ do + prop "prop_lookup" prop_lookup + prop "prop_fromToList" prop_fromToList + prop "prop_insertionPreserved" prop_insertionPreserved + prop "prop_nearest" prop_nearest + prop "prop_naiveWithinRange" prop_naiveWithinRange + prop "prop_fmap1" prop_fmap1 + prop "prop_fmap2" prop_fmap2 + prop "prop_fmap3" prop_fmap3 + prop "prop_depth_empty" prop_depth_empty + prop "prop_depth_upper_bound" prop_depth_upper_bound + prop "prop_size" prop_size diff --git a/tests/PropTests/OctreeTests/Exposed.hs b/tests/PropTests/OctreeTests/Exposed.hs new file mode 100644 --- /dev/null +++ b/tests/PropTests/OctreeTests/Exposed.hs @@ -0,0 +1,94 @@ +{- | + Module : PropTests.OctreeTests.Exposed + Copyright : Copyright (c) 2016 Michal J. Gajda + License : BSD3 + + Maintainer : Michael Litchard + Stability : experimental + Portability: not portable + + This module provides the property tests for exposed octree functions. + +-} + + +{-# LANGUAGE ScopedTypeVariables #-} +module PropTests.OctreeTests.Exposed + ( prop_lookup + , prop_fromToList + , prop_insertionPreserved + , prop_nearest + , prop_naiveWithinRange + , prop_fmap1 + , prop_fmap2 + , prop_fmap3 + , prop_depth_empty + , prop_depth_upper_bound + , prop_size + ) where + +import Prelude hiding(lookup) +import Data.List(sort, sortBy) + +import Data.Vector.Class +import Control.Arrow(second) + +import Data.Octree.Internal +import Data.Octree() -- test that interface module is not broken +import PropTests.Common + +prop_lookup :: [(Vector3, Int)] -> Bool +prop_lookup l = all isIn l + where + ot = fromList l + isIn x = lookup ot (fst x) == Just x + +prop_fromToList :: [(Vector3, Int)] -> Bool +prop_fromToList l = sort l == (sort . toList . fromList $ l) + +prop_insertionPreserved :: [(Vector3, Int)] -> Bool +prop_insertionPreserved l = + sort l == (sort . toList . foldr insert (Leaf []) $ l) + +prop_nearest :: [(Vector3, Int)] -> Vector3 -> Bool +prop_nearest l pt = nearest (fromList l) pt == naiveNearest pt l + +prop_naiveWithinRange r l pt = + naiveWithinRange r pt l == testPoints + where + testPoints = + sort . map fst . (\o -> withinRange o r pt) . fromList . tuplify pt $ l + +tuplify pt = map (\a -> (a, dist pt a)) + +compareDistance pt (a,_) (b,_) = compare (dist pt a) (dist pt b) + +naiveNearest pt [] = Nothing +naiveNearest pt l = Just $ head byDist + where byDist = sortBy (compareDistance pt) l + +naiveWithinRange r pt = sort . filter withinRange + where withinRange p = dist pt p <= r + +-- unfortunately there is no Arbitrary for (a -> b) +-- since generic properties are quite common, I wonder how to force Quickcheck to default something reasonable? +prop_fmap1,prop_fmap2 :: [(Vector3, Int)] -> Bool +prop_fmap1 = genericPropertyFmap (+1) +prop_fmap2 = genericPropertyFmap (*2) +prop_fmap3 = genericPropertyFmap (show :: Int -> String) + +genericPropertyFmap f l = + (sort . map (Control.Arrow.second f) $ l) == fmapTest + where fmapTest = sort . toList . fmap f . fromList $ l + +prop_depth_empty = depth (Leaf []) == 0 + +prop_depth_upper_bound :: [(Vector3, Int)] -> Bool +prop_depth_upper_bound l = + depth ot <= max 0 (ceiling . logBase 2 . realToFrac $ size) -- worst splitting ratio possible when we take midpoint (and inputs are colinear) + where + ot = fromList l + size = length l + +prop_size :: [(Vector3, Int)] -> Bool +prop_size l = size (fromList l) == length l diff --git a/tests/PropTests/OctreeTests/Internal.hs b/tests/PropTests/OctreeTests/Internal.hs new file mode 100644 --- /dev/null +++ b/tests/PropTests/OctreeTests/Internal.hs @@ -0,0 +1,143 @@ +{- | + Module : PropTests.OctreeTests.Internal + Copyright : Copyright (c) 2016 Michal J. Gajda + License : BSD3 + + Maintainer : Michael Litchard + Stability : experimental + Portability: not portable + + This module provides tests for internal helper functions: + +-} + +{-# LANGUAGE ScopedTypeVariables #-} +module PropTests.OctreeTests.Internal + ( prop_depth + , prop_cmp1 + , prop_cmp2 + , prop_stepDescription + , prop_octantDistanceNoGreaterThanInterpointDistance0 + , prop_octantDistanceNoGreaterThanInterpointDistance + , prop_octantDistanceNoGreaterThanInterpointDistanceZero + , prop_octantDistanceNoGreaterThanCentroidDistance + , prop_pickClosest + ) where + +import Prelude hiding(lookup) +import Data.List(sort, sortBy) + +import Control.Arrow(second) + +import Data.Octree.Internal +import Data.Octree() -- test that interface module is not broken + +import PropTests.Common + +-- for easier testing +origin :: Vector3 +origin = 0 + +prop_depth [] = True -- fudging , what is true fix? +prop_depth a = + (depth oct <= ((+1) . ceiling $ expectedDepth)) && + (depth oct >= ((\a -> a-1) . floor $ expectedDepth)) + where + expectedDepth = + (logBase 8 :: Double -> Double) . fromIntegral . length $ a + oct :: Octree Int = fromList a + +prop_cmp1 a b = cmp a b == joinStep (dx >= 0, dy >= 0, dz >= 0) + where Vector3 dx dy dz = a - b + +prop_cmp2 a = cmp a origin == joinStep (dx >= 0, dy >= 0, dz >= 0) + where Vector3 dx dy dz = a + +prop_stepDescription a b = + splitStep (cmp a b) == (v3x a >= v3x b, v3y a >= v3y b, v3z a >= v3z b) + +prop_octantDistanceNoGreaterThanInterpointDistance0 ptA ptB = triangleInequality + where + triangleInequality = octantDistance' aptA (cmp ptB origin) <= dist aptA ptB + aptA = abs ptA + +prop_octantDistanceNoGreaterThanInterpointDistance ptA ptB vp = + triangleInequality + where + triangleInequality = octantDistance (ptA - vp) (cmp ptB vp) <= dist ptA ptB + sameOctant = cmp ptA vp == cmp ptB vp + +prop_octantDistanceNoGreaterThanInterpointDistanceZero ptA ptB = + triangleInequality + where + triangleInequality = octantDistance ptA (cmp ptB origin) <= dist ptA ptB + sameOctant = cmp ptA origin == cmp ptB origin + +prop_octantDistanceNoGreaterThanInterpointDistanceZero0 ptA ptB = + triangleInequality + where + triangleInequality = octantDistance aptA (cmp ptB origin) <= dist aptA ptB + sameOctant = cmp aptA origin == cmp ptB origin + aptA = abs ptA + +prop_octantDistanceNoGreaterThanCentroidDistance pt vp = all testFun allOctants + where testFun odir = octantDistance (pt - vp) odir <= dist pt vp + +prop_splitByPrime splitPt pt = + (unLeaf . octreeStep ot . cmp pt $ splitPt) == [arg] + where + ot = splitBy' Leaf splitPt [arg] + arg = (pt, dist pt splitPt) + + +prop_pickClosest :: (Eq a) => [(Vector3, a)] -> Vector3 -> Bool +prop_pickClosest l pt = pickClosest pt l == naiveNearest pt l + +-- | These are tests for exposed functions: + +prop_lookup l = all isIn l + where + ot = fromList l + isIn x = lookup ot (fst x) == Just x + +prop_fromToList l = sort l == (sort . toList . fromList $ l) +prop_insertionPreserved l = + sort l == (sort . toList . foldr insert (Leaf []) $ l) +prop_nearest l pt = nearest (fromList l) pt == naiveNearest pt l +prop_naiveWithinRange r l pt = + naiveWithinRange r pt l == testPoints + where + testPoints = + sort . map fst . (\o -> withinRange o r pt) . fromList . tuplify pt $ l + +tuplify pt = map (\a -> (a, dist pt a)) + +compareDistance pt (a,_) (b,_) = compare (dist pt a) (dist pt b) + +naiveNearest pt [] = Nothing +naiveNearest pt l = Just $ head byDist + where byDist = sortBy (compareDistance pt) l + +naiveWithinRange r pt = sort . filter withinRange + where withinRange p = dist pt p <= r + +-- unfortunately there is no Arbitrary for (a -> b) +-- since generic properties are quite common, I wonder how to force Quickcheck to default something reasonable? +prop_fmap1,prop_fmap2 :: [(Vector3, Int)] -> Bool +prop_fmap1 = genericPropertyFmap (+1) +prop_fmap2 = genericPropertyFmap (*2) +prop_fmap3 = genericPropertyFmap (show :: Int -> String) + +genericPropertyFmap f l = + (sort . map (Control.Arrow.second f) $ l) == testFmap + where testFmap = sort . toList . fmap f . fromList $ l + +prop_depth_empty = depth (Leaf []) == 0 + +prop_depth_upper_bound l = + depth ot <= max 0 (ceiling . logBase 2 . realToFrac $ size) -- worst splitting ratio possible when we take midpoint (and inputs are colinear) + where + ot = fromList l + size = length l + +prop_size l = size (fromList l) == length l