packages feed

computational-geometry (empty) → 0.1.0

raw patch · 17 files changed

+1343/−0 lines, 17 filesdep +ansi-wl-pprintdep +basedep +containerssetup-changedbinary-added

Dependencies added: ansi-wl-pprint, base, containers, lens-family-core, linear, protolude, vector

Files

+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) 2017, Maksymilian Owsianny++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 Maksymilian Owsianny nor the names of other+      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.
+ README.md view
@@ -0,0 +1,16 @@+[![Linux Build Status][linux-build-icon]][linux-build]++# Computational Geometry++Collection of algorithms in Computational Geometry, specifically in the context+of procedural graphics generation. This is very much a work in progress.++Currently I'm working on set operations of polytopes. You can read more about+that in [This Blog Post][blog-post].++![Set Operations Example][setops3d]++[linux-build-icon]: https://img.shields.io/travis/MaxOw/computational-geometry/master.svg?label=build+[linux-build]: https://travis-ci.org/MaxOw/computational-geometry+[blog-post]: https://MaxOw.github.io/posts/computational-geometry-set-operations-on-polytopes.html+[setops3d]: images/setops3d.gif
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ computational-geometry.cabal view
@@ -0,0 +1,57 @@+Name          : computational-geometry+Version       : 0.1.0+Synopsis      : Collection of algorithms in Computational Geometry.+License       : BSD3+License-File  : LICENSE+Author        : Maksymilian Owsianny+Maintainer    : Maksymilian.Owsianny@gmail.com+Bug-Reports   : https//github.com/MaxOw/computational-geometry/issues+Category      : Graphics, Math+Build-Type    : Simple+Cabal-Version : >= 1.18++Description:+  Collection of algorithms in Computational Geometry.++Extra-Source-Files : README.md, images/*.png, images/*.gif+Extra-Doc-Files    : images/*.png, images/*.gif++Source-Repository head+  type:     git+  location: https://github.com/MaxOw/computational-geometry.git++Library+  default-language : Haskell2010+  hs-source-dirs   : src+  ghc-options      : -O2 -Wall -Wincomplete-uni-patterns++  exposed-modules:+    Data.EqZero+    Geometry.Plane.General+    Geometry.SetOperations+    Geometry.SetOperations.Types+    Geometry.SetOperations.Volume+    Geometry.SetOperations.Merge+    Geometry.SetOperations.BSP+    Geometry.SetOperations.CrossPoint+    Geometry.SetOperations.Facet+    Geometry.SetOperations.Clip+    Geometry.SetOperations.BRep++  default-extensions:+    NoImplicitPrelude+    DoAndIfThenElse+    LambdaCase+    MultiWayIf+    TupleSections+    OverloadedStrings++  build-depends : base >= 4.5 && < 5.0+                , protolude+                , containers+                , vector+                , linear++                , lens-family-core+                , ansi-wl-pprint+
+ images/set-operation-examples.png view

binary file changed (absent → 16660 bytes)

+ images/setops3d.gif view

binary file changed (absent → 354575 bytes)

+ src/Data/EqZero.hs view
@@ -0,0 +1,41 @@+{-# Language MultiParamTypeClasses #-}+{-# Language FlexibleInstances #-}+--------------------------------------------------------------------------------+-- |+-- Module     : Data.EqZero+-- Copyright  : (C) 2017 Maksymilian Owsianny+-- License    : BSD-style (see LICENSE)+-- Maintainer : Maksymilian.Owsianny@gmail.com+--+--------------------------------------------------------------------------------+module Data.EqZero where++import Protolude+import Linear.Epsilon (nearZero)+import Foreign.C (CFloat, CDouble)++--------------------------------------------------------------------------------++-- | Convenient universal zero equality predicate that warps to zero within some+-- epsilon for floating point numbers.+class EqZero a where+    eqZero :: a -> Bool++{-+instance (Num a, Eq a) => EqZero Exact a where+    eqZero _ = (==0)++instance Epsilon a => EqZero NonExact a where+    eqZero _ = nearZero+-}++-- | Greater or equal to zero predicate.+geqZero :: (EqZero n, Ord n, Num n) => n -> Bool+geqZero n = eqZero n || n > 0++instance EqZero Float    where eqZero = nearZero+instance EqZero Double   where eqZero = nearZero+instance EqZero CFloat   where eqZero = nearZero+instance EqZero CDouble  where eqZero = nearZero+instance EqZero Rational where eqZero = (==0)+
+ src/Geometry/Plane/General.hs view
@@ -0,0 +1,158 @@+{-# Language MultiParamTypeClasses #-}+{-# Language FlexibleInstances #-}+--------------------------------------------------------------------------------+-- |+-- Module     : Geometry.Plane.General+-- Copyright  : (C) 2017 Maksymilian Owsianny+-- License    : BSD-style (see LICENSE)+-- Maintainer : Maksymilian.Owsianny@gmail.com+--+-- General representation of a plane. Plane in the General Form is Hession+-- Normal Form scaled by an arbitrary non-zero scalar.+--+--------------------------------------------------------------------------------+module Geometry.Plane.General+    ( Plane (..)+    , Plane2, Plane3+    , Plane2D, Plane3D++    , MakePlane (..)+    , unsafeMakePlane+    , flipPlane++    , collinear+ -- , coincidence, coorientation++    , PlanesRelation (..), Incidence (..), Orientation (..)+    , planesRelation+    , isParallel++    ) where++import Protolude hiding (zipWith, zero)+import Data.Maybe (fromJust)+import qualified Data.List as List+import Linear+-- import Linear.Solve+import Linear.Affine (Point, (.-.))+import qualified Linear.Affine as Point+import Data.EqZero++-- | Internally Plane is represented as a pair (sN, sO) where N is a normal+-- vector of a plane O is the distance of that plane from the origin and s is an+-- arbitrary non-zero scalar.+data Plane v n = Plane+   { planeVector :: !(v n)+   , planeLast   :: !n+   } deriving (Eq, Ord, Show)++type Plane2 = Plane V2+type Plane3 = Plane V3++type Plane2D = Plane V2 Double+type Plane3D = Plane V3 Double++instance (NFData (v n), NFData n) => NFData (Plane v n) where+    rnf (Plane vs l) = rnf vs `seq` rnf l++-- | Flip plane orientation.+flipPlane :: (Functor v, Num n) => Plane v n -> Plane v n+flipPlane (Plane v n) = Plane (fmap negate v) (negate n)++class MakePlane v n where+    -- | Make plane from vector of points. Returns Nothing if vectors between+    -- points are linearly dependent+    makePlane :: v (Point v n) -> Maybe (Plane v n)++instance (Num n, Eq n) => MakePlane V3 n where+    makePlane (V3 p1 p2 p3)+        | n == zero = Nothing+        | otherwise = Just $ Plane n d+        where+        n = cross (p2 .-. p1) (p3 .-. p1)+        d = negate $ dot n $ unPoint p1++-- | Assumes that points form a valid plane (i.e. vectors between all points are+-- linearly independent).+unsafeMakePlane :: MakePlane v n => v (Point v n) -> Plane v n+unsafeMakePlane = fromJust . makePlane++{-+makePlane :: (Applicative v, Solve v n, Num n)+    => v (Point v n) -> Maybe (Plane v n)+-- makePlane ps = Plane <$> solve ups (pure 1) <*> pure 1+makePlane ps = uncurry Plane <$> solve ups (pure 1)+    where+    ups = fmap unPoint ps++-- | Assumes that points form a valid plane (i.e. vectors between all points are+-- linearly independent).+unsafeMakePlane :: (Applicative v, Solve v n, Num n)+    => v (Point v n) -> Plane v n+-- unsafeMakePlane ps = Plane (fromJust $ solve ups (pure 1)) 1+-- unsafeMakePlane ps = Plane v d+unsafeMakePlane ps = case solve ups (pure 1) of+    Just (v, d) -> Plane v d+    Nothing     -> error "Bla" -- . toS $ List.unlines $ map show ps+    where+    -- Just (v, d) = solve ups (pure 1)+    ups = fmap unPoint ps+-}++-- | Convert point to a vector.+unPoint :: Point v n -> v n+unPoint (Point.P x) = x++--------------------------------------------------------------------------------++-- | Test whether two vectors are collinear.+collinear :: (Foldable v, Num n, EqZero n) => v n -> v n -> Bool+collinear v w = all f $ combinations 2 $ zipWith (,) v w+    where+    f [(a, b), (c, d)] = eqZero $ a*d - b*c+    f _                = False -- To silence exhaustiveness checker++-- | All n-combinations of a given list.+combinations :: Int -> [a] -> [[a]]+combinations k is+    | k <= 0    = [ [] ]+    | otherwise = [ x:r | x:xs <- tails is, r <- combinations (k-1) xs ]++-- | Zip two `Foldable` structures to a list with a given function.+zipWith :: Foldable f => (a -> b -> c) -> f a -> f b -> [c]+zipWith f a b = List.zipWith f (toList a) (toList b)++-- | Test co-incidence of two planes assuming collinearity.+coincidence :: (Foldable v, Num n, EqZero n) => Plane v n -> Plane v n -> Bool+coincidence (Plane v1 d1) (Plane v2 d2) = all f $ zipWith (,) v1 v2+    where+    f (x1, x2) = eqZero $ x1*d2 - x2*d1++-- | Test co-orientation of two assuming collinearity.+coorientation :: (Foldable v, Num n, Ord n, EqZero n)+    => Plane v n -> Plane v n -> Bool+coorientation (Plane v1 d1) (Plane v2 d2)+    = all geqZero $ d1*d2 : zipWith (*) v1 v2++--------------------------------------------------------------------------------++data PlanesRelation = Parallel Incidence Orientation | Crossing deriving Show+data Incidence      = CoIncident |  NonIncident                 deriving Show+data Orientation    = CoOriented | AntiOriented                 deriving Show++-- | Relate two planes on Parallelism, Incidence and Orientation.+planesRelation :: (Foldable v, Num n, Ord n, EqZero n)+    => Plane v n -> Plane v n -> PlanesRelation+planesRelation p1@(Plane v1 _) p2@(Plane v2 _)+    | collinear v1 v2 = Parallel incidence orientation+    | otherwise       = Crossing+    where+    incidence   = bool  NonIncident CoIncident $ coincidence   p1 p2+    orientation = bool AntiOriented CoOriented $ coorientation p1 p2++isParallel :: (Foldable v, Num n, Ord n, EqZero n)+    => Plane v n -> Plane v n -> Bool+isParallel a b = case planesRelation a b of+    Parallel _ _ -> True+    Crossing     -> False+
+ src/Geometry/SetOperations.hs view
@@ -0,0 +1,167 @@+{-# Language ConstraintKinds #-}+{-# OPTIONS_HADDOCK prune #-}+--------------------------------------------------------------------------------+-- |+-- Module     : Geometry.SetOperations+-- Copyright  : (C) 2017 Maksymilian Owsianny+-- License    : BSD-style (see LICENSE)+-- Maintainer : Maksymilian.Owsianny@gmail.com+--+-- Set Operations of Polytopes. You can read about implementation details of+-- this algorithm in a dedicated <MaxOw.github.io/posts/computational-geometry-set-operations-on-polytopes.html Blog Post>.+--+-- Small example:+--+-- > test :: SetOperation -> Double -> PolyT3D+-- > test op t = fromVolume $ merge op boxA boxB+-- >     where+-- >     boxA = cube+-- >     boxB = toVolume $ Poly3 (papply tr <$> ps) is+-- >     Poly3 ps is = cubePoly3+-- >     tr = translation (V3 (sin (t*0.3) * 0.3) 0.2 0.3)+-- >        <> aboutX (t*20 @@ deg)+-- >        <> aboutY (t*3  @@ deg)+--+-- Rendered:+--+-- <<images/setops3d.gif>>+--+--------------------------------------------------------------------------------+module Geometry.SetOperations (++    -- * Base Functionality+      Volume, emptyVolume+    , toVolume, fromVolume++    , SetOperation (..)+    , merge, merges++    -- * Selected Merge Operations+    , union, unions+    , intersection, intersections+    , difference, differences++    -- * Conversion from/to BReps+    , FromPolytopeRep+    , ToPolytopeRep++    , Poly3 (..)+    , PolyT3 (..)++    -- * Primitives+    , cubePoly3, cube++    -- * Specializations/Synonyms+    , toVolume3D+    , fromVolume3D+    , Volume2D, Volume3D++    , Poly3D+    , PolyT3D+    , Merge++    ) where++import Protolude++import Linear+import Linear.Affine (Point)+import qualified Linear.Affine as Point++import qualified Data.Vector as T++import Geometry.SetOperations.Types+import Geometry.SetOperations.Volume+import Geometry.SetOperations.Clip+import Geometry.SetOperations.BRep++--------------------------------------------------------------------------------++-- | Convert an arbitrary polytope boundary representation into a Volume.+toVolume :: (FromPolytopeRep p b v n, Clip b v n, Functor v, Num n)+    => p v n -> Volume b v n+toVolume = makeVolume . fromPolytopeRep++-- | Recover a boundary representation of a Volume.+fromVolume :: ToPolytopeRep p b v n => Volume b v n -> p v n+fromVolume = toPolytopeRep . volumeFacets++-- | Convert a simple 3-BRep polyhedron to a Volume.+toVolume3D :: Poly3D -> Volume3D+toVolume3D = toVolume++-- | Reconstruct a triangulated 3-BRep from a Volume.+fromVolume3D :: Volume3D -> PolyT3D+fromVolume3D = fromVolume++--------------------------------------------------------------------------------++type Merge b v n = (Clip b v n, Functor v, Num n)++-- | Merge two Volumes under a specified Set Operation.+merge :: Merge b v n+      => SetOperation -> Volume b v n -> Volume b v n -> Volume b v n+merge = mergeVolumes++-- | Merges list of Volumes under a specified Set Operation. Empty list equals+-- empty set.+merges :: Merge b v n => SetOperation -> [Volume b v n] -> Volume b v n+merges _  []     = emptyVolume+merges op (v:vs) = foldl' (merge op) v vs+-- As to not leak memory on folding just a strict left fold is not enough. The+-- merge operation also needs to be strict since it operates on record with lazy+-- fields of spine lazy structures. Should I change representation to strict or+-- just deepseq it here? TODO: Fix this.++--------------------------------------------------------------------------------++-- | Union of two volumes. Convenience synonym for `merge Union`+union :: Merge b v n => Volume b v n -> Volume b v n -> Volume b v n+union = merge Union++-- | Union of list of volumes.+unions :: Merge b v n => [Volume b v n] -> Volume b v n+unions = merges Union++-- | Intersection of two volumes.+intersection :: Merge b v n => Volume b v n -> Volume b v n -> Volume b v n+intersection = merge Intersection++-- | Intersection of list of volumes.+intersections :: Merge b v n => [Volume b v n] -> Volume b v n+intersections = merges Intersection++-- | Difference between two volumes.+difference :: Merge b v n => Volume b v n -> Volume b v n -> Volume b v n+difference = merge Difference++-- | Subtract list of volumes from a given volume.+differences :: Merge b v n => Volume b v n -> [Volume b v n] -> Volume b v n+differences = foldl' (merge Difference)++--------------------------------------------------------------------------------++p3 :: a -> a -> a -> Point V3 a+p3 x y z = Point.P $ V3 x y z++-- | Cube represented as a denormalized list of polygons.+cubePoly3 :: Poly3D+cubePoly3 = Poly3 (T.fromList ps) is+    where+    ps = map (subtract 0.5) $+       [ p3 0 0 0, p3 1 0 0, p3 1 0 1, p3 0 0 1+       , p3 0 1 0, p3 1 1 0, p3 1 1 1, p3 0 1 1 ]++    is = map reverse+         [ [ 0, 1, 2, 3 ]+         , [ 1, 5, 6, 2 ]+         , [ 3, 2, 6, 7 ]+         , [ 0, 3, 7, 4 ]+         , [ 7, 6, 5, 4 ]+         , [ 0, 4, 5, 1 ]+         ]++-- | Cube volume.+cube :: Volume3D+cube = toVolume cubePoly3+
+ src/Geometry/SetOperations/BRep.hs view
@@ -0,0 +1,115 @@+{-# Language MultiParamTypeClasses #-}+{-# Language TypeSynonymInstances #-}+{-# Language FlexibleInstances #-}+--------------------------------------------------------------------------------+-- |+-- Module     : Geometry.SetOperations.BRep+-- Copyright  : (C) 2017 Maksymilian Owsianny+-- License    : BSD-style (see LICENSE)+-- Maintainer : Maksymilian.Owsianny@gmail.com+--+-- Boundary representations for conversion to and from BSP/Volumes.+--+--------------------------------------------------------------------------------+module Geometry.SetOperations.BRep+    ( FromPolytopeRep (..)+    , ToPolytopeRep   (..)++    , Poly3 (..), Poly3D+    , PolyT3 (..), PolyT3D+    ) where++import Protolude+import Linear.Affine (Point)+import Linear+import qualified Data.Map as Map++import Data.EqZero++-- import qualified Data.Vector.Generic as Vector+import Data.Vector.Generic ((!))+import qualified Data.Vector as T++import Geometry.Plane.General+import Geometry.SetOperations.Facet+import Geometry.SetOperations.CrossPoint+import Geometry.SetOperations.Clip++-- | Convert from polytope to a list of Facets.+class FromPolytopeRep p b v n where+    fromPolytopeRep :: p v n -> [Facet b v n]++-- | Convert from list of Facets to a polytope boundary representation.+class ToPolytopeRep p b v n where+    toPolytopeRep :: [Facet b v n] -> p v n++--------------------------------------------------------------------------------++-- | Indexed 3-BRep as a list of convex polygons. Continent as a way to+-- introduce new base shapes into the constructive geometry context.+data Poly3 v n = Poly3 (T.Vector (Point v n)) [[Int]]+type Poly3D = Poly3 V3 Double++instance ( MakePlane v n, Eq (v n), Foldable v, Applicative v, R3 v+         , Num n, Ord n, EqZero n+         ) => FromPolytopeRep Poly3 (FB3 v n) v n where+    fromPolytopeRep = makeFacets3++{-# SPECIALIZE makeFacets3 :: Poly3D -> [Facet3D] #-}++-- I assume valid indexes for now, without checks.+-- Will need to make it safe in the future.+-- There is also assumption that each point is shared by 3 planes+-- and that each eadge is shared by 2 planes.+makeFacets3 :: (MakePlane v n, Foldable v, Applicative v, R3 v, Ord n, EqZero n)+    => (Num n, Eq (v n))+    => Poly3 v n -> [Facet (FB3 v n) v n]+makeFacets3 (Poly3 ps is) = zipWith Facet planes boundries+    where+    points = map (map (ps!)) is+    planes = map (\(a:b:c:_) -> unsafeMakePlane $ vec3 a b c) points++    mkPlaneEdge (p, es) = map (,[p]) es++    edges    = map (map mkOrdPair . edges2) is+    edgesMap = Map.fromListWith (<>) $ concatMap mkPlaneEdge $ zip planes edges++    edgePlanePairs = map (mapMaybe (flip Map.lookup edgesMap)) edges+    edgePlanes     = zipWith edgeOnly planes edgePlanePairs+    edgeOnly p es  = map (\(a:b:_) -> if p == a then b else a) es++    uniqueCrossPoints = fmap toCrossPoint ps+    crossPoints       = map (map (uniqueCrossPoints!)) is++    boundries = zipWith (\a b -> zip a b) crossPoints edgePlanes++data OrdPair a = OrdPair !a !a deriving (Show, Eq, Ord)+mkOrdPair :: Ord a => (a, a) -> OrdPair a+mkOrdPair (a, b) = if a > b then OrdPair a b else OrdPair b a++{-# INLINE edges2 #-}+edges2 :: [a] -> [(a,a)]+edges2 as = zip as (drop 1 $ cycle as)++--------------------------------------------------------------------------------++-- | Simple direct 3-BRep as a list of triangles. Useful as an output after+-- performing specified set operations of the base shapes for rendering.+newtype PolyT3 v n = PolyT3 [ [Point v n] ]++type PolyT3D = PolyT3 V3 Double++instance ToPolytopeRep PolyT3 (FB3 v n) v n where+    toPolytopeRep fs = PolyT3 (concatMap f fs)+      where+      f (Facet _ bd) = tris $ map (getPoint . fst) bd++tris :: [a] -> [[a]]+tris ps = take triNum $ concat $ zipWith mkTri pps rps+    where+    triNum = length ps - 2+    pps    = egs ps+    rps    = egs $ reverse ps+    egs xs = zip xs $ drop 1 xs+    mkTri (a,b) (n,m) = [[a, m, n], [m, a, b]]+
+ src/Geometry/SetOperations/BSP.hs view
@@ -0,0 +1,130 @@+{-# Language PatternSynonyms   #-}+{-# Language DeriveFunctor     #-}+{-# Language OverloadedStrings #-}+--------------------------------------------------------------------------------+-- |+-- Module     : Geometry.SetOperations.BSP+-- Copyright  : (C) 2017 Maksymilian Owsianny+-- License    : BSD-style (see LICENSE)+-- Maintainer : Maksymilian.Owsianny@gmail.com+--+--------------------------------------------------------------------------------+module Geometry.SetOperations.BSP+    ( BinaryTree (..)+    , LeafColor  (..)+    , swapColor++    , BSP+    , cmp+    , pattern In+    , pattern Out++    , constructBSP+    , splitWith+    , destructBinaryTree++    , prettyBSP, renderH, denormalizeBSP+    ) where++import Prelude (id)+import Protolude hiding ((<>))+import Data.Monoid ((<>))++import Lens.Family (over)+import Lens.Family.Stock (both)+-- import Control.Lens (over, both)++import Data.List (unzip)+import Data.IntMap (IntMap)+import qualified Data.IntMap as IntMap+import qualified Data.Map as Map++import Text.PrettyPrint.ANSI.Leijen hiding ((<>), (<$>), dot, empty)++-- import Geometry.Plane.General+import Geometry.SetOperations.Facet+import Geometry.SetOperations.Clip++--------------------------------------------------------------------------------++-- | Binary Tree parametrized by leafs and nodes+data BinaryTree l n+   = Node (BinaryTree l n) !n (BinaryTree l n)+   | Leaf !l+   deriving (Eq, Show, Functor)++instance Bifunctor BinaryTree where+    bimap f _ (Leaf x)     = Leaf (f x)+    bimap f g (Node l n r) = Node (bimap f g l) (g n) (bimap f g r)++data LeafColor = Green | Red deriving (Eq, Show)++{-# INLINE swapColor #-}+swapColor :: LeafColor -> LeafColor+swapColor Green = Red+swapColor Red   = Green++type BSP = BinaryTree LeafColor++-- | Complementary set+cmp :: BSP a -> BSP a+cmp = first swapColor++pattern In :: BSP a+pattern In  = Leaf Green++pattern Out :: BSP a+pattern Out = Leaf Red++--------------------------------------------------------------------------------++constructBSP :: Clip b v n => (Facet b v n -> c) -> [Facet b v n] -> BSP c+constructBSP _ []                     = Out+constructBSP f (facet@(Facet s _):fs) = case splitWith (splitFacet s) fs of+    ([], rs) -> Node In                  c (constructBSP f rs)+    (ls, []) -> Node (constructBSP f ls) c Out+    (ls, rs) -> Node (constructBSP f ls) c (constructBSP f rs)+    where+    c = f facet++splitWith :: (a -> (Maybe a, Maybe a)) -> [a] -> ([a], [a])+splitWith f = over both catMaybes . unzip . map f++destructBinaryTree :: BinaryTree l n -> [n]+destructBinaryTree = flip go []+    where+    go (Node l p r) = (p:) . go l . go r+    go _            = identity++--------------------------------------------------------------------------------+-- Pretty Printing - for debugging+--------------------------------------------------------------------------------++type Context k = k -> Doc++-- | Pretty print BSP tree to stdout.+prettyBSP :: (Ord f) => BSP f -> IO ()+prettyBSP bsp = putDoc $ renderH id int bspId <+> linebreak+    where+    (bspId, _) = denormalizeBSP bsp++-- | Render BSP into a horizontal tree with a given context.+renderH :: (Doc -> Doc) -> Context k -> BSP k -> Doc+renderH _ _ In  = dullcyan "✔"+renderH _ _ Out = red      "✗"+renderH ind k (Node left pivot right) = vcat+    [ dullblue (k pivot)+    , ind $ "├ " <> renderH (ind . ("│ "<>)) k left+    , ind $ "└ " <> renderH (ind . ("  "<>)) k right+    ]++-- | Denormalize BSP with integers at nodes and IntMap of values.+denormalizeBSP :: Ord n => BSP n -> (BSP Int, IntMap n)+denormalizeBSP bsp = (fmap f bsp, fsMap)+    where+    fs    = ordNub $ destructBinaryTree bsp+    isMap = Map.fromList $ zip fs [0..]+    fsMap = IntMap.fromList $ zip [0..] fs++    f p = Map.findWithDefault (-1) p isMap+
+ src/Geometry/SetOperations/Clip.hs view
@@ -0,0 +1,148 @@+{-# Language MultiParamTypeClasses #-}+{-# Language TypeSynonymInstances #-}+{-# Language FlexibleInstances #-}+{-# Language DefaultSignatures #-}+--------------------------------------------------------------------------------+-- |+-- Module     : Geometry.SetOperations.Clip+-- Copyright  : (C) 2017 Maksymilian Owsianny+-- License    : BSD-style (see LICENSE)+-- Maintainer : Maksymilian.Owsianny@gmail.com+--+--------------------------------------------------------------------------------+module Geometry.SetOperations.Clip+    ( Clip (..)+    , vec3+    ) where++import Data.Function (id)+import Data.List (zipWith3, unzip)+import Protolude+import Linear+import Lens.Family ((.~), over)+import Lens.Family.Stock (both)+-- import Control.Lens ((.~), over, both)++import Data.EqZero+import Geometry.Plane.General+import Geometry.SetOperations.Facet+import Geometry.SetOperations.CrossPoint++--------------------------------------------------------------------------------++class Clip b v n where+    clipFacet  :: Plane v n   -- ^ Clipping plane+               -> Facet b v n -- ^ Facet to clip+               -> Maybe (Facet b v n)++    splitFacet :: Plane v n   -- ^ Splitting plane+               -> Facet b v n -- ^ Facet to split+               -> (Maybe (Facet b v n), Maybe (Facet b v n))++    clipFacet  p f = fst $ splitFacet p f+    default splitFacet :: (Functor v, Num n)+                       => Plane v n -> Facet b v n+                       -> (Maybe (Facet b v n), Maybe (Facet b v n))+    splitFacet p f = (clipFacet p f, clipFacet (flipPlane p) f)++    {-# MINIMAL (clipFacet | splitFacet) #-}++--------------------------------------------------------------------------------++splitCoincident :: (Foldable v, Num n, Ord n, EqZero n)+    => Plane v n -> Facet b v n+    -> (Maybe (Facet b v n), Maybe (Facet b v n))+    -> (Maybe (Facet b v n), Maybe (Facet b v n))+splitCoincident h f@(Facet s _) othercase = case planesRelation h s of+    Parallel CoIncident   CoOriented -> (Just f, Nothing)+    Parallel CoIncident AntiOriented -> (Nothing, Just f)+    _ -> othercase++vec2 :: (R2 v, Applicative v) => n -> n -> v n+vec2 x y = pure x & _xy .~ (V2 x y)++instance+    ( MakeCrossPoint v n, R2 v, Applicative v+    , Foldable v, Num n, Ord n, EqZero n )+    => Clip (FB2 v n) v n where+    splitFacet h f@(Facet s (a, b)) = splitCoincident h f othercase+      where+      mc     = makeCrossPoint $ vec2 h s+      go x y = Just $ Facet s (x, y)++      othercase = table (orientation a h) (orientation b h)+      table P M = (mc >>= \c -> go a c, mc >>= \c -> go c b)+      table M P = (mc >>= \c -> go c b, mc >>= \c -> go a c)+      table P _ = (Just f, Nothing)+      table _ P = (Just f, Nothing)+      table M _ = (Nothing, Just f)+      table _ M = (Nothing, Just f)+      -- This last case is not needed and is only here for completeness.+      -- It could happen if someone wrongly created a facet with edge+      -- points not lying on the facet plane (line). In such case, that+      -- facet is simply discarded by the splitting function.+      table Z Z = (Nothing, Nothing)++--------------------------------------------------------------------------------++vec3 :: (R3 v, Applicative v) => n -> n -> n -> v n+vec3 x y z = pure x & _xyz .~ (V3 x y z)++instance+    ( MakeCrossPoint v n, R3 v, Applicative v+    , Foldable v, Num n, Ord n, EqZero n )+    => Clip (FB3 v n) v n where+    splitFacet h f@(Facet s ps) = splitCoincident h f othercase+      where+      mc v = makeCrossPoint $ vec3 s h v+      go ops@(_:_:_:_) = Just $ Facet s ops+      go _             = Nothing+      ss = map (flip orientation h . fst) ps++      othercase = over both go $ splitFast mc h ss ps++splitFast+    :: (p -> Maybe c)       -- ^ Make CrossPoint from V+    -> p                    -- ^ Clipping plane H+    -> [Sign]               -- ^ Points signs relative to H+    -> [(c, p)]             -- ^ Cross Boundry+    -> ([(c, p)], [(c, p)]) -- ^ Result+splitFast mkP h ss pvs+    | all (/= M) ss = (pvs, [])+    | all (/= P) ss = ([], pvs)+    | otherwise     = (compose outPlus, compose outMinus)+    where+    (outPlus, outMinus) = unzip $ zipWith3 table pvs ss (dropCycle 1 ss)++    table (p, v) P M = case mkP v of+        Nothing -> (mk1 (p, v), id)+        Just c  -> (mk2 (p, v) (c, h), mk1 (c, v))+    table (p, v) M P = case mkP v of+        Nothing -> (id, mk1 (p, v))+        Just c  -> (mk1 (c, v), mk2 (p, v) (c, h))++    table (p, v) Z M = (mk1 (p, v), mk1 (p, h))+    table (p, v) Z P = (mk1 (p, h), mk1 (p, v))++    table pv     P _ = (mk1 pv, id)+    table pv     M _ = (id, mk1 pv)++    table _      _ _ = (id, id) -- This case should never happen+    -- If it happens it means that it's a concave boundry.++{-# INLINE compose #-}+compose :: [([a] -> [a])] -> [a]+compose fs = foldr (.) id fs []++{-# INLINE mk1 #-}+mk1 :: a -> ([a] -> [a])+mk1 a   = (a:)++{-# INLINE mk2 #-}+mk2 :: a -> a -> ([a] -> [a])+mk2 a b = (a:) . (b:)++{-# INLINE dropCycle #-}+dropCycle :: Int -> [a] -> [a]+dropCycle n = drop n . cycle+
+ src/Geometry/SetOperations/CrossPoint.hs view
@@ -0,0 +1,99 @@+{-# Language MultiParamTypeClasses #-}+{-# Language FlexibleInstances     #-}+--------------------------------------------------------------------------------+-- |+-- Module     : Geometry.SetOperations.CrossPoint+-- Copyright  : (C) 2017 Maksymilian Owsianny+-- License    : BSD-style (see LICENSE)+-- Maintainer : Maksymilian.Owsianny@gmail.com+--+--------------------------------------------------------------------------------+module Geometry.SetOperations.CrossPoint+    ( Sign (..)+    , toSign++    , CrossPoint     (..)+    , MakeCrossPoint (..)+    , toCrossPoint++    ) where++import Protolude+import Linear+import Linear.Affine (Point)+import qualified Linear.Affine as Point++import Data.EqZero+import Geometry.Plane.General++--------------------------------------------------------------------------------++data Sign = M | Z | P deriving (Show, Eq)++toSign :: (EqZero n, Ord n, Num n) => n -> Sign+toSign x+    | eqZero x  = Z+    | x < 0     = M+    | otherwise = P++data CrossPoint v n = CP+   { orientation :: Plane v n -> Sign+   , getPoint    :: Point v n+   }++-- | Convert a point to CrossPoint+toCrossPoint :: (EqZero n, Foldable v, Num n, Ord n)+    => Point v n -> CrossPoint v n+toCrossPoint pt = CP orient pt+    where+    orient p = toSign . ((planeLast p) +) . sum+             $ zipWith (*) (toList $ planeVector p) (toList pt)++class MakeCrossPoint v n where+    makeCrossPoint :: v (Plane v n) -> Maybe (CrossPoint v n)++instance (Fractional n, Ord n, EqZero n) => MakeCrossPoint V2 n where+    makeCrossPoint planes+        | eqZero d2 = Nothing+        | otherwise = Just $ CP orient solved+        where+        V2     (Plane (V2 a b) c)+               (Plane (V2 d e) f) = planes+     -- orient (Plane (V2 g h) i) = toSign $ d2*(g*d0 - h*d1 + i*d2)+        orient (Plane (V2 g h) i) = toSign $ g*dd0 + h*dd1 + i++        dd0 = d2*d0+        dd1 = d2*d1++        d0 = b*f - c*e+        d1 = a*f - c*d+        d2 = a*e - b*d++        dd = 1/d2+        solved = Point.P $ V2 (dd*d0) (dd*d1)++instance (Fractional n, Ord n, EqZero n) => MakeCrossPoint V3 n where+    makeCrossPoint planes+        | eqZero d3 = Nothing+        | otherwise = Just $ CP orient solved+        where+        V3     (Plane (V3 a b c) d)+               (Plane (V3 e f g) h)+               (Plane (V3 i j k) l) = planes+        orient (Plane (V3 m n o) p) = toSign $ -d3*(m*d0 - n*d1 + o*d2 - p*d3)++        d0 = k*m1 - j*m0 + l*m2+        d1 = k*m3 - i*m0 + l*m4+        d2 = j*m3 - i*m1 + l*m5+        d3 = i*m2 - j*m4 + k*m5++        m0 = c*h - d*g+        m1 = b*h - d*f+        m2 = c*f - b*g+        m3 = a*h - d*e+        m4 = c*e - a*g+        m5 = b*e - a*f++        dd = 1/d3+        solved = Point.P $ V3 (-dd*d0) (dd*d1) (-dd*d2)+
+ src/Geometry/SetOperations/Facet.hs view
@@ -0,0 +1,39 @@+--------------------------------------------------------------------------------+-- |+-- Module     : Geometry.SetOperations.Facet+-- Copyright  : (C) 2017 Maksymilian Owsianny+-- License    : BSD-style (see LICENSE)+-- Maintainer : Maksymilian.Owsianny@gmail.com+--+--------------------------------------------------------------------------------+module Geometry.SetOperations.Facet+    ( Facet (..)+    , Facet2D, Facet3D+    , flipFacet+    , FB2, FB3++    ) where++import Protolude+import Linear (V2, V3)++import Geometry.Plane.General+import Geometry.SetOperations.CrossPoint++--------------------------------------------------------------------------------++data Facet b v n = Facet+   { facetPlane    :: Plane v n+   , facetBoundary :: b+   }++-- | Flip orientation of a facet.+flipFacet :: (Functor v, Num n) => Facet b v n -> Facet b v n+flipFacet (Facet p b) = Facet (flipPlane p) b++type FB3 v n = [(CrossPoint v n, Plane v n)]+type FB2 v n = (CrossPoint v n, CrossPoint v n)++type Facet2D = Facet (FB2 V2 Double) V2 Double+type Facet3D = Facet (FB3 V3 Double) V3 Double+
+ src/Geometry/SetOperations/Merge.hs view
@@ -0,0 +1,233 @@+{-# Language MultiParamTypeClasses       #-}+{-# Language TypeSynonymInstances        #-}+{-# Language FlexibleInstances           #-}+{-# OPTIONS_GHC -Wno-unused-matches      #-}+{-# OPTIONS_GHC -Wno-incomplete-patterns #-}+--------------------------------------------------------------------------------+-- |+-- Module     : Geometry.SetOperations.Merge+-- Copyright  : (C) 2017 Maksymilian Owsianny+-- License    : BSD-style (see LICENSE)+-- Maintainer : Maksymilian.Owsianny@gmail.com+--+-- Set Operations of Polytopes by BSP Merging.+--+--------------------------------------------------------------------------------+module Geometry.SetOperations.Merge+    ( BSP+    , BSP3D, BSP2D++    , Universe (..)+    , universePlanes, universeBox+    , splitRegion++    , mergeBSPs+    , trim++    , makeBSP+    , toBoundary+    ) where++import Protolude+import Prelude (id)++import Lens.Family (over)+import Lens.Family.Stock (both, _2)+-- import Control.Lens (over, both, _2)++import Data.Maybe (fromMaybe, fromJust)+import Linear++import Geometry.SetOperations.Types+import Geometry.SetOperations.BSP+import Geometry.SetOperations.Facet+import Geometry.SetOperations.CrossPoint+import Geometry.SetOperations.Clip+import Geometry.Plane.General+import Data.EqZero++type BSP2D = BSP Facet2D+type BSP3D = BSP Facet3D++--------------------------------------------------------------------------------++-- Arbitrary selected as sufficient by independent comity (not really).+universeSize :: Num n => n+universeSize = 500++clipPlanes :: Clip b v n => Facet b v n -> [Plane v n] -> Facet b v n+clipPlanes = foldr (\p f -> fromMaybe f $ clipFacet p f)++class Clip b v n => Universe b v n where+    -- | Turn plane into a Facet by clipping it by the universe box.+    makeFacet :: Plane v n -> Facet b v n++instance (Ord n, Fractional n, EqZero n) => Universe (FB2 V2 n) V2 n where+    makeFacet p = clipPlanes baseFacet ps+        where+        baseFacet = Facet p (a, b)+        Just a = makeCrossPoint (V2 p pa)+        Just b = makeCrossPoint (V2 p pb)+        (pa:pb:ps) = filter (not . isParallel p) universePlanes++instance (Ord n, Fractional n, EqZero n) => Universe (FB3 V3 n) V3 n where+    makeFacet p = Facet p es+        where+        ps = filter (not . isParallel p) universePlanes+        es = zipWith mkBd ps $ drop 1 $ cycle ps+        mkBd a b = (fromJust . makeCrossPoint $ V3 p a b, b)++-- | Planes bounding the UniverseBox.+universePlanes :: (Applicative v, Traversable v, Num n) => [Plane v n]+universePlanes = positive ++ negative+    where+    toPlane v    = Plane v universeSize+    positive     = map toPlane (basisFor $ pure 0)+    negative     = map flipPlane positive++-- | List of facets bounding the Universe.+universeBox :: (Universe b v n, Applicative v, Traversable v, Num n)+    => [Facet b v n]+universeBox = map makeFacet universePlanes++-- | Split a region within a Universe bounded by a list of Facets.+splitRegion :: (Universe b v n, Functor v, Num n)+    => Plane v n -> [Facet b v n] -> ([Facet b v n], [Facet b v n])+splitRegion h fs = (flipFacet lid : plusC, lid : minusC)+    where+    (plusC, minusC) = splitWith (splitFacet h) fs+    lid = clipPlanes (makeFacet h) (map facetPlane fs)++{-+type Merge b v n =+    (Universe b v n, Applicative v, Traversable v, Num n, Ord n, EqZero n)+-}++-- | Perform a given SetOperation of two BSPs by merging+mergeBSPs+    :: (Universe b v n, Applicative v, Traversable v, Num n, Ord n, EqZero n)+    => SetOperation+    -> BSP (Facet b v n)+    -> BSP (Facet b v n)+    -> BSP (Facet b v n)+mergeBSPs op (Node treeL p treeR) nodeR@(Node _ f _) =+    collapse $ Node mTreeL p mTreeR+    where+    ff = facetPlane f+    pp = facetPlane p+    regions = splitRegion ff universeBox+    (partL, partR) = partitionBSP regions pp nodeR+    mTreeL = mergeBSPs op treeL partL+    mTreeR = mergeBSPs op treeR partR+mergeBSPs op s1 s2 = setOperation op s1 s2++partitionBSP+    :: (Universe b v n, Functor v, Foldable v, Num n, Ord n, EqZero n)+    => ([Facet b v n], [Facet b v n])+    -> Plane v n+    -> BSP (Facet b v n)+    -> (BSP (Facet b v n), BSP (Facet b v n))+partitionBSP _       _ (Leaf c)             = (Leaf c, Leaf c)+partitionBSP regions p (Node treeP f treeM) = case planesRelation p ff of+    Parallel CoIncident   CoOriented -> (treeP, treeM)+    Parallel CoIncident AntiOriented -> (treeM, treeP)+    othercase -> if+      | null regionPR -> (Node treeP f treeML, treeMR)+      | null regionMR -> (Node treePL f treeM, treePR)+      | null regionPL -> (treeML, Node treeP f treeMR)+      | null regionML -> (treePL, Node treePR f treeM)++      | otherwise     -> (Node treePL f treeML, Node treePR f treeMR)+    where+    ff = facetPlane f+    (treePL, treePR) = partitionBSP (regionPL, regionPR) p treeP+    (treeML, treeMR) = partitionBSP (regionML, regionMR) p treeM++    (regionP , regionM ) = regions+    (regionPL, regionPR) = splitRegion p regionP+    (regionML, regionMR) = splitRegion p regionM++setOperation :: SetOperation -> BSP a -> BSP a -> BSP a++setOperation    Union                  In   set = In+setOperation    Union                  Out  set = set+setOperation    Union                  set  In  = In+setOperation    Union                  set  Out = set++setOperation    Intersection           In   set = set+setOperation    Intersection           Out  set = Out+setOperation    Intersection           set  In  = set+setOperation    Intersection           set  Out = Out++setOperation    Difference             In   set = cmp set+setOperation    Difference             Out  set = Out+setOperation    Difference             set  In  = Out+setOperation    Difference             set  Out = set++setOperation    SymmetricDifference    In   set = cmp set+setOperation    SymmetricDifference    Out  set = set+setOperation    SymmetricDifference    set  In  = cmp set+setOperation    SymmetricDifference    set  Out = set++collapse :: BSP n -> BSP n+collapse (Node In  _ In ) = In+collapse (Node Out _ Out) = Out+collapse other            = other++isBoundary :: Clip b v n => BSP (Facet b v n) -> Facet b v n -> Bool+isBoundary In  _ = True+isBoundary Out _ = False+isBoundary (Node l s r) f = lcnd || rcnd+    where+    (lh, rh) = splitFacet (facetPlane s) f+    lcnd = fromMaybe False (isBoundary l <$> lh)+    rcnd = fromMaybe False (isBoundary r <$> rh)++-- | Optimize a resulting BSP after merging by removing superficial splitting+-- planes.+trim :: Clip b v n => BSP (Facet b v n) -> BSP (Facet b v n)+trim (Node Out f r)+    | isBoundary r f = Node Out f (trim r)+    | otherwise     = trim r+trim (Node l f Out)+    | isBoundary l f = Node (trim l) f Out+    | otherwise     = trim l+trim other = other++--------------------------------------------------------------------------------++-- | Make a BSP from a list of bounding facets.+makeBSP :: Clip b v n => [Facet b v n] -> BSP (Facet b v n)+makeBSP = constructBSP id++-- | Reconstruct boundary facets from the BSP.+toBoundary :: (Clip b v n, Functor v, Num n)+    => BSP (Facet b v n) -> [Facet b v n]+toBoundary bsp+    = removeColors+    . map (over _2 flipFacet)+    . applyColors+    $ destructBinaryTree bsp+    where+    applyColors xs = go xs bsp []+        where+        go [] _   = id+        go fs In  = foldr (\f cs -> ((True , f):) . cs) id fs+        go fs Out = foldr (\f cs -> ((False, f):) . cs) id fs+        go fs (Node l s r) = go ls l . go rs r+            where+            sp = facetPlane s+            (ls, rs) = splitWith (splitFacet sp) fs++    removeColors xs = go xs bsp []+        where+        go [] _   = id+        go fs In  = foldr (\(a,b) cs -> if not a then (b:) . cs else cs) id fs+        go fs Out = foldr (\(a,b) cs -> if     a then (b:) . cs else cs) id fs+        go fs (Node l s r) = go ls l . go rs r+            where+            (ls, rs) = splitWith coloredSplit fs+            sp = facetPlane s+            coloredSplit (b, f) = over both (fmap (b,)) $ splitFacet sp f++
+ src/Geometry/SetOperations/Types.hs view
@@ -0,0 +1,24 @@+--------------------------------------------------------------------------------+-- |+-- Module     : Geometry.SetOperations.Types+-- Copyright  : (C) 2017 Maksymilian Owsianny+-- License    : BSD-style (see LICENSE)+-- Maintainer : Maksymilian.Owsianny@gmail.com+--+--------------------------------------------------------------------------------+module Geometry.SetOperations.Types+    ( SetOperation (..)+    ) where++--------------------------------------------------------------------------------++-- | Four basic set operations:+--+-- <<images/set-operation-examples.png>>++data SetOperation+   = Union+   | Intersection+   | Difference+   | SymmetricDifference+
+ src/Geometry/SetOperations/Volume.hs view
@@ -0,0 +1,84 @@+--------------------------------------------------------------------------------+-- |+-- Module     : Geometry.SetOperations.Volume+-- Copyright  : (C) 2017 Maksymilian Owsianny+-- License    : BSD-style (see LICENSE)+-- Maintainer : Maksymilian.Owsianny@gmail.com+--+-- Set Operations of Polytopes by Boundary Filtering.+--+--------------------------------------------------------------------------------+module Geometry.SetOperations.Volume+    ( Volume (..)+    , makeVolume+    , emptyVolume+    , mergeVolumes++    , Volume2D, Volume3D+    ) where++import Protolude+import Linear (V2, V3)++import Geometry.SetOperations.Merge+import Geometry.SetOperations.Types+import Geometry.SetOperations.BSP+import Geometry.SetOperations.Facet+import Geometry.SetOperations.Clip+import Geometry.Plane.General++--------------------------------------------------------------------------------++-- | Volume, currently represented as a list of Facets and a BSP Tree.+data Volume b v n = Volume+   { volumeFacets :: [Facet b v n]+   , volumeTree   :: BSP (Plane v n)+   }++type Volume2D = Volume (FB2 V2 Double) V2 Double+type Volume3D = Volume (FB3 V3 Double) V3 Double++-- | Construct Volume from a list of Facets representing it's boundary.+makeVolume :: Clip b v n => [Facet b v n] -> Volume b v n+makeVolume fs = Volume fs (constructBSP facetPlane fs)++-- | Empty volume.+emptyVolume :: Volume b v n+emptyVolume = Volume [] Out++--------------------------------------------------------------------------------++{-# SPECIALIZE+  mergeVolumes :: SetOperation -> Volume2D -> Volume2D -> Volume2D #-}+{-# SPECIALIZE+  mergeVolumes :: SetOperation -> Volume3D -> Volume3D -> Volume3D #-}++-- | Merge two Volumes under a specified Set Operation.+mergeVolumes :: (Clip b v n, Functor v, Num n)+    => SetOperation -> Volume b v n -> Volume b v n -> Volume b v n+mergeVolumes op volumeA volumeB = case op of+    Difference          -> filterBoth isOut    isInFlip+    Intersection        -> filterBoth isIn     isIn+    Union               -> filterBoth isOut    isOut+    SymmetricDifference -> filterBoth isEither isEither+    where+    isInFlip x fs = case x of Red -> []; Green -> map flipFacet fs+    isIn     x fs = case x of Red -> []; Green -> fs+    isOut    x fs = case x of Red -> fs; Green -> []+    isEither x fs = case x of Red -> fs; Green -> map flipFacet fs++    Volume facetsA treeA = volumeA+    Volume facetsB treeB = volumeB++    filterBoth f g = makeVolume $+        filterWith f facetsA treeB <>+        filterWith g facetsB treeA++    filterWith _ [] _ = []+    filterWith f fs t = case t of+        Leaf x             -> f x fs+        Node treeL p treeR ->+            filterWith f partL treeL <>+            filterWith f partR treeR+            where (partL, partR) = splitWith (splitFacet p) fs+