diff --git a/Data/Octree.hs b/Data/Octree.hs --- a/Data/Octree.hs +++ b/Data/Octree.hs @@ -1,5 +1,5 @@ module Data.Octree(Octree, - Vector3(..), + Vector3, dist, fromList, toList, lookup, @@ -14,4 +14,3 @@ import Prelude hiding(lookup) import Data.Octree.Internal - diff --git a/Data/Octree/Internal.hs b/Data/Octree/Internal.hs --- a/Data/Octree/Internal.hs +++ b/Data/Octree/Internal.hs @@ -1,6 +1,7 @@ {-# LANGUAGE ScopedTypeVariables, DisambiguateRecordFields #-} {-# LANGUAGE DeriveFunctor, DeriveFoldable, DeriveTraversable #-} -module Data.Octree.Internal(Vector3(..), dist, +{-# LANGUAGE DeriveGeneric, DeriveAnyClass #-} +module Data.Octree.Internal(Vector3(..), dist, v3x, v3y, v3z, Octree(..), lookup, nearest, withinRange, fromList, toList, insert, delete, deleteBy, -- internal ODir(..), @@ -15,9 +16,14 @@ import Data.Foldable (Foldable (foldr)) import Data.Traversable(Traversable(..)) -import Data.Vector.V3 -import Data.Vector.Class +import GHC.Generics(Generic) +import Control.DeepSeq(NFData) +import Control.Lens.Getter +import Linear.Metric +import Linear.V3 +import Linear + --import Text.Show import Prelude hiding(lookup, foldr) import qualified Data.List as List (delete, deleteBy, sortBy) @@ -27,24 +33,28 @@ import Test.QuickCheck.All(quickCheckAll) import Test.QuickCheck.Arbitrary --- | norm of a vector -norm :: Vector3 -> Double -norm a = sqrt (a `vdot` a) +type Vector3 = V3 Double +v3x, v3y, v3z :: V3 Double -> Double +v3x v = v ^. _x +v3y v = v ^. _y +v3z v = v ^. _z + -- | distance between two vectors -dist :: Vector3 -> Vector3 -> Double -dist u v = norm (u - v) +dist :: V3 Double -> V3 Double -> Double +dist = distance -- | Datatype for nodes within Octree. -data Octree a = Node { split :: Vector3, +data Octree a = Node { split :: V3 Double, nwu, nwd, neu, ned, swu, swd, seu, sed :: Octree a } | - Leaf { unLeaf :: [(Vector3, a)] } deriving (Show, Functor, Foldable, Traversable) + Leaf { unLeaf :: [(V3 Double, a)] } + deriving (Show, Functor, Foldable, Traversable, Generic, NFData) -- | 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) -- | Internal method that gives octant of a first vector relative to the second vector as a center. -cmp :: Vector3 -> Vector3 -> ODir +cmp :: V3 Double -> V3 Double -> ODir cmp ca cb = joinStep (cx, cy, cz) where cx = v3x ca >= v3x cb cy = v3y ca >= v3y cb @@ -78,7 +88,7 @@ -- 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' :: V3 Double -> ODir -> Double octantDistance' dp NEU = 0.0 -- adjacent by plane octantDistance' dp NWU = v3x dp @@ -101,12 +111,12 @@ -- | Finds a minimum bounds for a distance between a given point -- | in relative coordinates and a given octant. -octantDistance :: Vector3 -> ODir -> Double +octantDistance :: V3 Double -> 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 :: V3 Double -> ODir -> ODir toggle dp odir = joinStep ((v3x dp >= 0) `xor` not u, (v3y dp >= 0) `xor` not v, @@ -114,20 +124,20 @@ 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 :: V3 Double -> [(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)] -> ([(Vector3, a)], - [(Vector3, a)], - [(Vector3, a)], - [(Vector3, a)], - [(Vector3, a)], - [(Vector3, a)], - [(Vector3, a)], - [(Vector3, a)]) +splitBy :: V3 Double -> [(V3 Double, a)] -> ([(V3 Double, a)], + [(V3 Double, a)], + [(V3 Double, a)], + [(V3 Double, a)], + [(V3 Double, a)], + [(V3 Double, a)], + [(V3 Double, a)], + [(V3 Double, a)]) splitBy _splitPoint [] = ([], [], [], [], [], [], [], []) splitBy splitPoint ((pt@(coord, a)):aList) = case i of @@ -157,16 +167,16 @@ leafLimit = 16 -- | Creates an Octree from a list of (index, payload) tuples. -fromList :: [(Vector3, a)] -> Octree a +fromList :: [(V3 Double, a)] -> Octree a fromList aList = if length aList <= leafLimit then Leaf aList - else let splitPoint :: Vector3 = massCenter aList + else let splitPoint :: V3 Double = 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' :: ([(V3 Double, a)] -> Octree a1)-> V3 Double-> [(V3 Double, a)]-> Octree a1 splitBy' f splitPoint aList = Node { split = splitPoint, nwu = tnwu, nwd = tnwd, @@ -182,7 +192,7 @@ -- | Internal method that prepends contents of the given subtree to a list -- | given as argument. -toList' :: Octree t -> [(Vector3, t)] -> [(Vector3, t)] +toList' :: Octree t -> [(V3 Double, t)] -> [(V3 Double, t)] toList' (Leaf l ) tmp = l ++ tmp toList' (Node { nwu = a, nwd = b, @@ -194,12 +204,12 @@ 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 :: Octree t -> [(V3 Double, 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 :: V3 Double -> Octree a -> [ODir] pathTo pt (Leaf _) = [] pathTo pt node = aStep : pathTo pt (octreeStep node aStep) where aStep = cmp pt (split node) @@ -220,7 +230,7 @@ -- | 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 :: (V3 Double, 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) @@ -228,7 +238,7 @@ -- | Deletes a point from an Octree. -- | NOTE: If there are duplicate points, it only deletes one of them. -delete :: (Eq a) => (Vector3, a) -> Octree a -> Octree a +delete :: (Eq a) => (V3 Double, a) -> Octree a -> Octree a delete tup@(pt, _) ot = applyByPath delete' path ot where path = pathTo pt ot delete' (Leaf l) = fromList (List.delete tup l) @@ -236,7 +246,7 @@ -- | Deletes a point from an Octree with the provided equality check. -- | NOTE: If there are duplicate points, it only deletes one of them. -deleteBy :: ((Vector3, a) -> (Vector3, a) -> Bool) -> (Vector3, a) -> Octree a -> Octree a +deleteBy :: ((V3 Double, a) -> (V3 Double, a) -> Bool) -> (V3 Double, a) -> Octree a -> Octree a deleteBy eq tup@(pt, _) ot = applyByPath deleteBy' path ot where path = pathTo pt ot deleteBy' (Leaf l) = fromList (List.deleteBy eq tup l) @@ -244,7 +254,7 @@ -- | 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' :: V3 Double -> Octree a -> [(ODir, Double, [(V3 Double, a)])] candidates' pt (Leaf l) = [] candidates' pt node = map findCandidates . List.sortBy compareDistance . octantDistances $ pt - split node where @@ -252,12 +262,12 @@ 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 :: Octree a -> V3 Double -> Maybe (V3 Double, 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 :: Octree a -> V3 Double -> Maybe (V3 Double, a) nearest (Leaf l) pt = pickClosest pt l nearest node pt = selectFrom candidates where candidates = map findCandidate . List.sortBy compareDistance . octantDistances $ pt - split node @@ -279,7 +289,7 @@ 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 :: V3 Double -> [(V3 Double, t)] -> Maybe (V3 Double, 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 @@ -287,7 +297,7 @@ else vb -- | Returns all points within Octree that are within a given distance from argument. -withinRange :: Octree a -> Scalar -> Vector3 -> [(Vector3, a)] +withinRange :: Octree a -> Double -> V3 Double -> [(V3 Double, 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 diff --git a/Octree.cabal b/Octree.cabal --- a/Octree.cabal +++ b/Octree.cabal @@ -1,5 +1,5 @@ name: Octree -version: 0.5.4.4 +version: 0.6.0.0 stability: beta homepage: https://github.com/mgajda/octree package-url: http://hackage.haskell.org/package/octree @@ -29,7 +29,11 @@ location: git@github.com:mgajda/octree.git Library - build-depends: base>=4.0 && < 4.11, AC-Vector >= 2.3.0, QuickCheck >= 2.4.0 + build-depends: base>=4.0 && < 4.11, + deepseq >= 1.4.1.0, + lens >= 4.13.0, + linear >= 1.20.0, + QuickCheck >= 2.4.0 exposed-modules: Data.Octree other-modules: Data.Octree.Internal exposed: True @@ -37,13 +41,19 @@ Test-suite test_Octree Type: exitcode-stdio-1.0 - Build-depends: base>=4.0 && < 4.11, AC-Vector >= 2.3.0, QuickCheck >= 2.8.0 + Build-depends: base>=4.0 && < 4.11, linear >= 1.20.0, QuickCheck >= 2.8.0, + deepseq >= 1.4.1.0, + lens >= 4.13.0 + other-modules: Data.Octree + Data.Octree.Internal Main-is: tests/test_Octree.hs ---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.11, AC-Vector >= 2.3.0, QuickCheck >= 2.4.0, markdown-unlit --- ghc-options: -pgmL markdown-unlit - +benchmark bench-octree + type: exitcode-stdio-1.0 + other-modules: Data.Octree.Internal + main-is: tests/benchmark.hs + build-depends: base >=4.0 && <4.11, deepseq, ghc-prim, + linear >= 1.20.0, QuickCheck >= 2.8.0, + lens >= 4.13.0, + criterion + extensions: ScopedTypeVariables diff --git a/changelog b/changelog --- a/changelog +++ b/changelog @@ -1,4 +1,8 @@ -*-Changelog-*- +0.6.0.0 Apr 2018 + * Switched 3D vectors to `linear`. + * Added first benchmark. + 0.5.4.4 Dec 2017 * Remove doctest since it is fragile (depends on markdown-unlit being on path.) diff --git a/tests/benchmark.hs b/tests/benchmark.hs new file mode 100644 --- /dev/null +++ b/tests/benchmark.hs @@ -0,0 +1,27 @@ +{-# LANGUAGE FlexibleInstances #-} +module Main(main) where + +import Criterion +import Criterion.Main +import Test.QuickCheck(generate, vectorOf, Arbitrary(..)) + +import Linear.V3 +import Data.Octree.Internal + +sizes = [100, 1000, 10000] + +instance {-# OVERLAPPING #-} Arbitrary (Vector3, Int) where + arbitrary = ((,) <$> (V3 <$> arbitrary <*> arbitrary <*> arbitrary) + <*> arbitrary) + +testData :: Int -> IO [(Vector3, Int)] +testData i = generate $ vectorOf i arbitrary + +main = defaultMain [ + bgroup "insert" [ + let i=10^p + in bench (show i) $ nfIO (foldr insert (Leaf []) + <$> testData i) + | p <- [1..4] + ] + ] diff --git a/tests/test_Octree.hs b/tests/test_Octree.hs --- a/tests/test_Octree.hs +++ b/tests/test_Octree.hs @@ -1,4 +1,5 @@ {-# LANGUAGE TemplateHaskell, ScopedTypeVariables #-} +{-# LANGUAGE FlexibleInstances #-} module Main where import Data.Octree.Internal @@ -12,23 +13,14 @@ import Test.QuickCheck.Gen (elements) import Test.QuickCheck.Modifiers (NonEmptyList (..)) -import Data.Vector.Class +import Linear import Control.Arrow(second) --- | 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 +instance Arbitrary (V3 Double) where arbitrary = do x <- arbitrary y <- arbitrary z <- arbitrary - return $ Vector3 x y z + return $ V3 x y z data ElementListWithElement a = ElementListWithElement [(Vector3, a)] (Vector3, a) deriving (Show) @@ -48,18 +40,19 @@ prop_depth a = (depth oct <= ((+1) . ceiling $ expectedDepth)) && (depth oct >= ((\a -> a-1) . floor $ expectedDepth)) where - expectedDepth = (logBase 8 :: Double -> Double) . fromIntegral . length $ a + expectedDepth = max 0 + $ (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 + where V3 dx dy dz = a - b prop_cmp2 a = cmp a origin == joinStep (dx >= 0, dy >= 0, dz >= 0) - where Vector3 dx dy dz = a + where V3 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 +prop_octantDistanceNoGreaterThanInterpointDistance0 ptA ptB = triangleInequality where triangleInequality = octantDistance' aptA (cmp ptB origin) <= dist aptA ptB aptA = abs ptA @@ -137,4 +130,3 @@ return [] main = $quickCheckAll -