diff --git a/bampy/src/slicer/split_surface.rs b/bampy/src/slicer/split_surface.rs
index 7407d0d..e04a9aa 100644
--- a/bampy/src/slicer/split_surface.rs
+++ b/bampy/src/slicer/split_surface.rs
@@ -7,42 +7,50 @@ pub fn split_surface(mut triangles: Vec<Triangle<f64>>) -> Vec<Mesh<f64>> {
     while let Some(triangle) = triangles.pop() {
         let mut surface = vec![triangle];
         let mut bvh = Bvh::build(&mut surface);
-        let mut aabb = surface[0].aabb;
-        triangles.retain_mut(|triangle| {
-            let mut stack = Vec::<usize>::new();
-            stack.push(0);
-            while let Some(i) = stack.pop() {
-                match bvh.nodes[i] {
-                    BvhNode::Node {
-                        parent_index: _,
-                        child_l_index,
-                        child_l_aabb,
-                        child_r_index,
-                        child_r_aabb,
-                    } => {
-                        if triangle.has_point_in_aabb(&child_l_aabb) {
-                            stack.push(child_l_index);
+        let mut aabb = surface[0].aabb.clone();
+
+        loop {
+            let mut mutated = false;
+            triangles.retain_mut(|triangle| {
+                let mut stack = Vec::<usize>::new();
+                stack.push(0);
+                while let Some(i) = stack.pop() {
+                    match bvh.nodes[i] {
+                        BvhNode::Node {
+                            parent_index: _,
+                            child_l_index,
+                            child_l_aabb,
+                            child_r_index,
+                            child_r_aabb,
+                        } => {
+                            if triangle.has_point_in_aabb(&child_l_aabb) {
+                                stack.push(child_l_index);
+                            }
+                            if triangle.has_point_in_aabb(&child_r_aabb) {
+                                stack.push(child_r_index);
+                            }
                         }
-                        if triangle.has_point_in_aabb(&child_r_aabb) {
-                            stack.push(child_r_index);
-                        }
-                    }
-                    BvhNode::Leaf {
-                        parent_index: _,
-                        shape_index,
-                    } => {
-                        if triangle.connected_with_triangle(surface[shape_index]) {
-                            surface.push(*triangle);
-                            let index = surface.len() - 1;
-                            bvh.add_shape(&mut surface, index);
-                            aabb.join_mut(&triangle.aabb);
-                            return false;
+                        BvhNode::Leaf {
+                            parent_index: _,
+                            shape_index,
+                        } => {
+                            if triangle.connected_with_triangle(surface[shape_index]) {
+                                mutated = true;
+                                surface.push(*triangle);
+                                let index = surface.len() - 1;
+                                bvh.add_shape(&mut surface, index);
+                                aabb.join_mut(&triangle.aabb);
+                                return false;
+                            }
                         }
                     }
                 }
+                true
+            });
+            if !mutated {
+                break;
             }
-            true
-        });
+        }
 
         surfaces.push(Mesh {
             triangles: surface,
diff --git a/bampy/src/slicer/triangle.rs b/bampy/src/slicer/triangle.rs
index 923c511..54f8ed9 100644
--- a/bampy/src/slicer/triangle.rs
+++ b/bampy/src/slicer/triangle.rs
@@ -18,16 +18,19 @@ pub struct Triangle<T: SimdPartialOrd + Scalar + Copy> {
     pub aabb: Aabb<T, 3>,
 }
 
-fn vec_inside_aabb<T: SimdPartialOrd + Scalar + Copy + Float>(
+#[inline(always)]
+fn vec_inside_aabb<T: SimdPartialOrd + Scalar + Copy + Float + RelativeEq + approx::AbsDiffEq>(
     vec: &Vector3<T>,
     aabb: &Aabb<T, 3>,
 ) -> bool {
-    vec.x >= aabb.min.x
-        && vec.y >= aabb.min.y
-        && vec.z >= aabb.min.z
-        && vec.x <= aabb.max.x
-        && vec.y <= aabb.max.y
-        && vec.z <= aabb.max.z
+    macro_rules! within {
+        ($axis:ident) => {
+            ((vec.$axis >= aabb.min.$axis && vec.$axis <= aabb.max.$axis)
+                || relative_eq!(vec.$axis, aabb.min.$axis)
+                || relative_eq!(vec.$axis, aabb.max.$axis))
+        };
+    }
+    within!(x) && within!(y) && within!(z)
 }
 
 impl<T> Triangle<T>