diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -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.
diff --git a/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -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
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/computational-geometry.cabal b/computational-geometry.cabal
new file mode 100644
--- /dev/null
+++ b/computational-geometry.cabal
@@ -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
+
diff --git a/images/set-operation-examples.png b/images/set-operation-examples.png
new file mode 100644
Binary files /dev/null and b/images/set-operation-examples.png differ
diff --git a/images/setops3d.gif b/images/setops3d.gif
new file mode 100644
Binary files /dev/null and b/images/setops3d.gif differ
diff --git a/src/Data/EqZero.hs b/src/Data/EqZero.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/EqZero.hs
@@ -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)
+
diff --git a/src/Geometry/Plane/General.hs b/src/Geometry/Plane/General.hs
new file mode 100644
--- /dev/null
+++ b/src/Geometry/Plane/General.hs
@@ -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
+
diff --git a/src/Geometry/SetOperations.hs b/src/Geometry/SetOperations.hs
new file mode 100644
--- /dev/null
+++ b/src/Geometry/SetOperations.hs
@@ -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
+
diff --git a/src/Geometry/SetOperations/BRep.hs b/src/Geometry/SetOperations/BRep.hs
new file mode 100644
--- /dev/null
+++ b/src/Geometry/SetOperations/BRep.hs
@@ -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]]
+
diff --git a/src/Geometry/SetOperations/BSP.hs b/src/Geometry/SetOperations/BSP.hs
new file mode 100644
--- /dev/null
+++ b/src/Geometry/SetOperations/BSP.hs
@@ -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
+
diff --git a/src/Geometry/SetOperations/Clip.hs b/src/Geometry/SetOperations/Clip.hs
new file mode 100644
--- /dev/null
+++ b/src/Geometry/SetOperations/Clip.hs
@@ -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
+
diff --git a/src/Geometry/SetOperations/CrossPoint.hs b/src/Geometry/SetOperations/CrossPoint.hs
new file mode 100644
--- /dev/null
+++ b/src/Geometry/SetOperations/CrossPoint.hs
@@ -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)
+
diff --git a/src/Geometry/SetOperations/Facet.hs b/src/Geometry/SetOperations/Facet.hs
new file mode 100644
--- /dev/null
+++ b/src/Geometry/SetOperations/Facet.hs
@@ -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
+
diff --git a/src/Geometry/SetOperations/Merge.hs b/src/Geometry/SetOperations/Merge.hs
new file mode 100644
--- /dev/null
+++ b/src/Geometry/SetOperations/Merge.hs
@@ -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
+
+
diff --git a/src/Geometry/SetOperations/Types.hs b/src/Geometry/SetOperations/Types.hs
new file mode 100644
--- /dev/null
+++ b/src/Geometry/SetOperations/Types.hs
@@ -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
+
diff --git a/src/Geometry/SetOperations/Volume.hs b/src/Geometry/SetOperations/Volume.hs
new file mode 100644
--- /dev/null
+++ b/src/Geometry/SetOperations/Volume.hs
@@ -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
+