diff --git a/Example/SceneKitViewController.swift b/Example/SceneKitViewController.swift
index 85fe9b6f..b5bc56a7 100644
--- a/Example/SceneKitViewController.swift
+++ b/Example/SceneKitViewController.swift
@@ -32,6 +32,10 @@ class SceneKitViewController: UIViewController {
 //        let sphere = Mesh.sphere(slices: 120, material: UIColor.blue)
 //        let mesh = cube.subtracting(sphere).makeWatertight()
 
+//        let cube = Mesh.cube().inverted()
+//        let cube2 = Mesh.cube().inverted().translated(by: Vector(0.5, 0.5, 0.5))
+//        let mesh = cube.union(cube2).inverted()
+
 //        print("Time:", CFAbsoluteTimeGetCurrent() - start)
 //        print("Polygons:", mesh.polygons.count)
 //        print("Triangles:", mesh.triangulate().polygons.count)
diff --git a/Sources/BSP.swift b/Sources/BSP.swift
index 08e072da..cbb0f9f7 100644
--- a/Sources/BSP.swift
+++ b/Sources/BSP.swift
@@ -46,10 +46,14 @@ extension BSP {
     }
 
     init(_ mesh: Mesh, _ isCancelled: CancellationHandler) {
+        self.init(mesh.polygons, isConvex: mesh.isKnownConvex, isCancelled)
+    }
+
+    init(_ polygons: [Polygon], isConvex: Bool, _ isCancelled: CancellationHandler) {
         self.nodes = [BSPNode]()
-        self.isConvex = mesh.isKnownConvex
+        self.isConvex = isConvex
         self.isInverted = false
-        initialize(mesh.polygons, isCancelled)
+        initialize(polygons, isCancelled)
     }
 
     func clip(
@@ -138,12 +142,13 @@ private extension BSP {
     }
 
     mutating func initialize(_ polygons: [Polygon], _ isCancelled: CancellationHandler) {
+        guard !polygons.isEmpty else {
+            return
+        }
+
         var rng = DeterministicRNG()
 
         guard isConvex else {
-            guard !polygons.isEmpty else {
-                return
-            }
             let startPlane = polygons[0].plane
             // Randomly shuffle polygons to reduce average number of splits
             let polygons = polygons.shuffled(using: &rng)
diff --git a/Sources/Mesh+CSG.swift b/Sources/Mesh+CSG.swift
index cf73e570..6c4ecdd4 100644
--- a/Sources/Mesh+CSG.swift
+++ b/Sources/Mesh+CSG.swift
@@ -54,9 +54,40 @@ public extension Mesh {
     /// - Returns: A new mesh representing the union of the input meshes.
     func union(_ mesh: Mesh, isCancelled: CancellationHandler = { false }) -> Mesh {
         let intersection = bounds.intersection(mesh.bounds)
-        if intersection.isEmpty {
+        let absp = getBSP(isCancelled), bbsp = mesh.getBSP(isCancelled)
+        switch (absp.isInverted, bbsp.isInverted) {
+        case (false, false):
+            if intersection.isEmpty {
+                return Mesh(
+                    unchecked: polygons + mesh.polygons,
+                    bounds: bounds.union(mesh.bounds),
+                    isConvex: false,
+                    isWatertight: watertightIfSet.flatMap { isWatertight in
+                        mesh.watertightIfSet.map { $0 && isWatertight }
+                    },
+                    submeshes: [self, mesh]
+                )
+            }
+            var lhs: [Polygon] = [], rhs: [Polygon] = []
+            inParallel({
+                var aout: [Polygon]? = []
+                let ap = BSP(mesh, isCancelled).clip(
+                    boundsTest(intersection, polygons, &aout),
+                    .greaterThan,
+                    isCancelled
+                )
+                lhs = aout! + ap
+            }, {
+                var bout: [Polygon]? = []
+                let bp = BSP(self, isCancelled).clip(
+                    boundsTest(intersection, mesh.polygons, &bout),
+                    .greaterThanEqual,
+                    isCancelled
+                )
+                rhs = bout! + bp
+            })
             return Mesh(
-                unchecked: polygons + mesh.polygons,
+                unchecked: lhs + rhs,
                 bounds: bounds.union(mesh.bounds),
                 isConvex: false,
                 isWatertight: watertightIfSet.flatMap { isWatertight in
@@ -66,32 +97,47 @@ public extension Mesh {
                     mesh.submeshesIfEmpty.map { _ in [self, mesh] }
                 }
             )
-        }
-        var lhs: [Polygon] = [], rhs: [Polygon] = []
-        inParallel({
-            var aout: [Polygon]? = []
-            let ap = BSP(mesh, isCancelled).clip(
-                boundsTest(intersection, polygons, &aout),
+        case (true, true):
+            if intersection.isEmpty {
+                return .empty
+            }
+            var out: [Polygon]?
+            let ap = bbsp.clip(
+                boundsTest(intersection, polygons, &out),
                 .greaterThan,
                 isCancelled
             )
-            lhs = aout! + ap
-        }, {
-            var bout: [Polygon]? = []
-            let bp = BSP(self, isCancelled).clip(
-                boundsTest(intersection, mesh.polygons, &bout),
+            let bp = absp.clip(
+                boundsTest(intersection, mesh.polygons, &out),
                 .greaterThanEqual,
                 isCancelled
             )
-            rhs = bout! + bp
-        })
-        return Mesh(
-            unchecked: lhs + rhs,
-            bounds: bounds.union(mesh.bounds),
-            isConvex: false,
-            isWatertight: nil,
-            submeshes: nil // TODO: can this be preserved?
-        )
+            return Mesh(
+                unchecked: ap + bp,
+                bounds: bounds.union(mesh.bounds),
+                isConvex: false,
+                isWatertight: nil,
+                submeshes: nil // TODO: can this be preserved?
+            )
+        default:
+            let ap = bbsp.clip(
+                polygons,
+                .greaterThan,
+                isCancelled
+            )
+            let bp = absp.clip(
+                mesh.polygons,
+                .greaterThanEqual,
+                isCancelled
+            )
+            return Mesh(
+                unchecked: ap + bp,
+                bounds: bounds.union(mesh.bounds),
+                isConvex: false,
+                isWatertight: nil,
+                submeshes: nil // TODO: can this be preserved?
+            )
+        }
     }
 
     /// Efficiently forms a union from multiple meshes.
@@ -210,7 +256,6 @@ public extension Mesh {
             let (bp2, bp1) = absp.split(bp, .greaterThan, .lessThan, isCancelled)
             rhs = bout! + bp2 + bp1.inverted()
         })
-
         return Mesh(
             unchecked: lhs + rhs,
             bounds: nil, // TODO: is there a way to efficiently preserve this?
@@ -337,7 +382,7 @@ public extension Mesh {
         }
         var aout: [Polygon]? = []
         let ap = boundsTest(bounds.intersection(mesh.bounds), polygons, &aout)
-        let bsp = BSP(mesh, isCancelled)
+        let bsp = mesh.getBSP(isCancelled)
         let (outside, inside) = bsp.split(ap, .greaterThan, .lessThanEqual, isCancelled)
         let material = mesh.polygons.first?.material
         return Mesh(
@@ -408,11 +453,16 @@ public extension Mesh {
 
     /// Clip mesh to the specified plane and optionally fill sheared faces with specified material.
     /// - Parameters
-    ///   - plane: The plane to clip the mesh to
-    ///   - fill: The material to fill the sheared face(s) with.
+    ///   - plane: The plane to clip the mesh against.
+    ///   - fill: Optional material to fill the sheared face(s) with.
+    ///   - isCancelled: Callback used to cancel the operation.
     ///
     /// > Note: Specifying nil for the fill material will leave the sheared face unfilled.
-    func clip(to plane: Plane, fill: Material? = nil) -> Mesh {
+    func clip(
+        to plane: Plane,
+        fill: Material? = nil,
+        isCancelled: CancellationHandler = { false }
+    ) -> Mesh {
         guard !polygons.isEmpty else {
             return self
         }
diff --git a/Sources/Mesh.swift b/Sources/Mesh.swift
index 0a343ade..3432cd7c 100644
--- a/Sources/Mesh.swift
+++ b/Sources/Mesh.swift
@@ -339,7 +339,7 @@ public extension Mesh {
     /// - Returns: `true` if the point lies inside the mesh, and `false` otherwise.
     func containsPoint(_ point: Vector) -> Bool {
         guard isKnownConvex else {
-            return BSP(self) { false }.containsPoint(point)
+            return storage.getBSP { false }.containsPoint(point)
         }
         if !bounds.containsPoint(point) {
             return false
@@ -386,11 +386,16 @@ extension Mesh {
 
     var boundsIfSet: Bounds? { storage.boundsIfSet }
     var watertightIfSet: Bool? { storage.watertightIfSet }
+    var bspIfSet: Bool? { storage.watertightIfSet }
     var isKnownConvex: Bool { storage.isConvex }
     /// Note: we don't expose submeshesIfSet because it's unsafe to reuse
     var submeshesIfEmpty: [Mesh]? {
         storage.submeshesIfSet.flatMap { $0.isEmpty ? [] : nil }
     }
+
+    func getBSP(_ isCancelled: CancellationHandler) -> BSP {
+        storage.getBSP(isCancelled)
+    }
 }
 
 private extension Mesh {
@@ -437,6 +442,18 @@ private extension Mesh {
             return watertightIfSet!
         }
 
+        private(set) var bspIfSet: BSP?
+        func getBSP(_ isCancelled: CancellationHandler) -> BSP {
+            var bsp = bspIfSet
+            if bsp == nil {
+                bsp = BSP(polygons, isConvex: isConvex, isCancelled)
+                if !isCancelled() {
+                    bspIfSet = bsp
+                }
+            }
+            return bsp!
+        }
+
         private(set) var submeshesIfSet: [Mesh]?
         var submeshes: [Mesh] {
             if submeshesIfSet == nil {