diff --git a/Kinetic_space_partition/examples/Kinetic_space_partition/include/Parameters.h b/Kinetic_space_partition/examples/Kinetic_space_partition/include/Parameters.h
index 3a89c0e64f0..53682e3ce9a 100644
--- a/Kinetic_space_partition/examples/Kinetic_space_partition/include/Parameters.h
+++ b/Kinetic_space_partition/examples/Kinetic_space_partition/include/Parameters.h
@@ -59,11 +59,9 @@ namespace KSR {
       debug(false),
       // shape detection / shape regularization
       k_neighbors(12),
-      min_region_size(0),
+      min_region_size(100),
       max_octree_node_size(40),
       max_octree_depth(3),
-      distance_threshold(0),
-      angle_threshold(10),
       // partition
       k_intersections(1),
       enlarge_bbox_ratio(FT(11) / FT(10)),
diff --git a/Kinetic_space_partition/include/CGAL/KSP/debug.h b/Kinetic_space_partition/include/CGAL/KSP/debug.h
index ee34ce72c6f..d2f8b408ea3 100644
--- a/Kinetic_space_partition/include/CGAL/KSP/debug.h
+++ b/Kinetic_space_partition/include/CGAL/KSP/debug.h
@@ -915,41 +915,14 @@ void dump_volumes(const DS& data, const std::string tag = std::string()) {
   }
 }
 
-/*
-template<typename KSP>
-void dump_volumes_ksp(const KSP& ksp, const std::string tag = std::string()) {
-  using Point_3 = typename KSP::Point_3;
-  using Index = typename KSP::Index;
-  std::vector< std::vector<Point_3> > polygons;
-  std::vector<Color> colors;
+template<typename K>
+void dump_polygon(const std::vector<typename K::Point_3>& pts, const std::string& filename) {
+  Saver<K> saver;
+  std::vector<std::vector<typename K::Point_3> > pts2;
+  pts2.push_back(pts);
 
-  std::vector<Index> faces_of_volume;
-  std::vector<Point_3> pts_of_face;
-
-  Saver<typename KSP::Kernel> saver;
-  for (std::size_t i = 0; i < ksp.number_of_volumes(); ++i) {
-    faces_of_volume.clear();
-    polygons.clear();
-    colors.clear();
-
-    const auto color = saver.get_idx_color(i);
-    ksp.faces(i, std::back_inserter(faces_of_volume));
-
-    colors.clear();
-    polygons.clear();
-    for (const Index& f : faces_of_volume) {
-      pts_of_face.clear();
-      ksp.vertices(f, std::back_inserter(pts_of_face));
-
-      polygons.push_back(pts_of_face);
-      colors.push_back(color);
-    }
-
-    const std::string file_name = tag + std::to_string(i);
-    saver.export_polygon_soup_3(polygons, colors, file_name);
-  }
+  saver.export_polygon_soup_3(pts2, filename);
 }
-*/
 
 void dump_polygon(const std::vector<CGAL::Epick::Point_3>& pts, const std::string& filename) {
   Saver<CGAL::Epick> saver;
diff --git a/Kinetic_space_partition/include/CGAL/KSP/utils.h b/Kinetic_space_partition/include/CGAL/KSP/utils.h
index 7a39034a8a8..d8e99ddddc2 100644
--- a/Kinetic_space_partition/include/CGAL/KSP/utils.h
+++ b/Kinetic_space_partition/include/CGAL/KSP/utils.h
@@ -68,7 +68,7 @@ decltype(auto) distance(const Point_d& p, const Point_d& q) {
   using Traits = typename Kernel_traits<Point_d>::Kernel;
   using FT = typename Traits::FT;
   const FT sq_dist = CGAL::squared_distance(p, q);
-  return static_cast<FT>(CGAL::sqrt(CGAL::to_double(sq_dist)));
+  return static_cast<FT>(CGAL::approximate_sqrt(sq_dist));
 }
 
 // Project 3D point onto 2D plane.
@@ -94,7 +94,7 @@ inline const Vector_d normalize(const Vector_d& v) {
   using FT = typename Traits::FT;
   const FT dot_product = CGAL::abs(v * v);
   //CGAL_assertion(dot_product != FT(0));
-  return v / static_cast<FT>(CGAL::sqrt(CGAL::to_double(dot_product)));
+  return v / static_cast<FT>(CGAL::approximate_sqrt(dot_product));
 }
 
 
diff --git a/Kinetic_space_partition/include/CGAL/KSP_2/Data_structure.h b/Kinetic_space_partition/include/CGAL/KSP_2/Data_structure.h
index f6a2c9f78d6..5ec82ee8cee 100644
--- a/Kinetic_space_partition/include/CGAL/KSP_2/Data_structure.h
+++ b/Kinetic_space_partition/include/CGAL/KSP_2/Data_structure.h
@@ -38,7 +38,7 @@ public:
   typedef typename Kernel::Line_2 Line_2;
   typedef typename Kernel::Segment_2 Segment_2;
 
-  typedef Support_line<Kernel> Support_line_DS;
+  typedef Support_line_2<Kernel> Support_line_DS;
 
   typedef CGAL::KSP_2::internal::Vertex<FT> Vertex;
   typedef CGAL::KSP_2::internal::Segment Segment;
diff --git a/Kinetic_space_partition/include/CGAL/KSP_3/Data_structure.h b/Kinetic_space_partition/include/CGAL/KSP_3/Data_structure.h
index f02d7ba90aa..5a64606321a 100644
--- a/Kinetic_space_partition/include/CGAL/KSP_3/Data_structure.h
+++ b/Kinetic_space_partition/include/CGAL/KSP_3/Data_structure.h
@@ -45,6 +45,7 @@ public:
   using Face_event = typename Support_plane::Face_event;
 
   using FT = typename Kernel::FT;
+  using IkFT = typename Intersection_kernel::FT;
   using Point_2 = typename Kernel::Point_2;
   using IkPoint_2 = typename Intersection_kernel::Point_2;
   using Point_3 = typename Kernel::Point_3;
@@ -54,6 +55,7 @@ public:
   using Segment_3 = typename Kernel::Segment_3;
   using IkSegment_3 = typename Intersection_kernel::Segment_3;
   using Vector_2 = typename Kernel::Vector_2;
+  using IkVector_2 = typename Intersection_kernel::Vector_2;
   using Direction_2 = typename Kernel::Direction_2;
   using IkDirection_2 = typename Intersection_kernel::Direction_2;
   using Triangle_2  = typename Kernel::Triangle_2;
@@ -240,7 +242,8 @@ private:
   Reconstructed_model m_reconstructed_model;
 
 public:
-  Data_structure(const Parameters& parameters, const std::string &prefix) : to_exact(), from_exact(), m_parameters(parameters), m_prefix(prefix) { }
+  Data_structure(const Parameters& parameters, const std::string &prefix) : to_exact(), from_exact(), m_parameters(parameters), m_prefix(prefix) {
+  }
 
   template<typename Type1, typename Type2, typename ResultType>
   static bool intersection(const Type1& t1, const Type2& t2, ResultType& result) {
@@ -365,31 +368,40 @@ public:
     m_support_planes.resize(number_of_items);
   }
 
-  FT calculate_edge_intersection_time(std::size_t sp_idx, IEdge edge, Face_event &event) {
+  IkFT calculate_edge_intersection_time(std::size_t sp_idx, IEdge edge, Face_event &event) {
     // Not need to calculate for border edges.
     if (m_intersection_graph.iedge_is_on_bbox(edge))
       return 0;
 
+    // Count faces
+    std::size_t numfaces = 0;
+    for (std::size_t i = 0; i < m_support_planes.size(); i++)
+      numfaces += m_support_planes[i].data().mesh.number_of_faces();
+
     Support_plane& sp = m_support_planes[sp_idx];
 
+    To_exact to_exact;
+
     Point_2 centroid = sp.data().centroid;
+    IkPoint_2 centroid2 = to_exact(sp.data().centroid);
 
     typename Intersection_graph::Kinetic_interval& kinetic_interval = m_intersection_graph.kinetic_interval(edge, sp_idx);
 
-/*
-    if (kinetic_interval.size() != 0) {
-      int a;
-      a = 3;
-    }*/
-
     Point_2 s = sp.to_2d(from_exact(point_3(m_intersection_graph.source(edge))));
+    IkPoint_2 s2 = sp.to_2d(point_3(m_intersection_graph.source(edge)));
     Point_2 t = sp.to_2d(from_exact(point_3(m_intersection_graph.target(edge))));
+    IkPoint_2 t2 = sp.to_2d(point_3(m_intersection_graph.target(edge)));
     Vector_2 segment = t - s;
-    FT segment_length = sqrt(segment * segment);
+    IkVector_2 segment2 = t2 - s2;
+    FT segment_length = CGAL::approximate_sqrt(segment * segment);
+    IkFT segment_length2 = CGAL::approximate_sqrt(segment2.squared_length());
     CGAL_assertion(segment_length > 0);
     segment = segment / segment_length;
+    segment2 = segment2 / segment_length2;
     Direction_2 to_source(s - centroid);
     Direction_2 to_target(t - centroid);
+    IkDirection_2 to_source2(s2 - centroid2);
+    IkDirection_2 to_target2(t2 - centroid2);
 
     const std::size_t uninitialized = static_cast<std::size_t>(-1);
     std::size_t source_idx = uninitialized;
@@ -419,8 +431,9 @@ public:
 
     std::size_t num;
 
-    Vector_2 tt = to_target.vector(), ts = to_source.vector();
-    bool ccw = (tt.x() * ts.y() - tt.y() * ts.x()) < 0;
+    // Vector_2 tt = to_target.vector(), ts = to_source.vector();
+    IkVector_2 tt2 = to_target2.vector(), ts2 = to_source2.vector();
+    bool ccw = (tt2.x() * ts2.y() - tt2.y() * ts2.x()) < 0;
 
     // Check whether the segment is cw or ccw oriented.
     if (!ccw) {
@@ -431,6 +444,10 @@ public:
       Point_2 tmp_p = s;
       s = t;
       t = tmp_p;
+
+      IkPoint_2 tmp_p2 = s2;
+      s2 = t2;
+      t2 = tmp_p2;
     }
 
     if (source_idx <= target_idx)
@@ -438,12 +455,13 @@ public:
     else
       num = (sp.data().original_directions.size() + target_idx - source_idx);
 
-    std::vector<FT> time(num);
-    std::vector<Point_2> intersections(num);
-    std::vector<FT> intersections_bary(num);
+    std::vector<IkFT> time(num);
+    std::vector<IkPoint_2> intersections(num);
+    std::vector<IkFT> intersections_bary(num);
 
     // Shooting rays to find intersection with line of IEdge
-    typename Intersection_kernel::Line_2 l = sp.to_2d(m_intersection_graph.line_3(edge));
+    typename Intersection_kernel::Line_3 l3 = m_intersection_graph.line_3(edge);
+    const typename Intersection_kernel::Line_2 l = sp.to_2d(l3);
     for (std::size_t i = 0; i < num; i++) {
       std::size_t idx = (i + source_idx) % sp.data().original_directions.size();
       const auto result = CGAL::intersection(l, sp.data().original_rays[idx]);
@@ -452,14 +470,17 @@ public:
         continue;
       }
       IkPoint_2 p;
-      if (CGAL::assign(p, result)) {
-        FT l = CGAL::sqrt(sp.data().original_vectors[idx].squared_length());
 
-        double l2 = CGAL::to_double((p - sp.data().original_rays[idx].point(0)).squared_length());
-        time[i] = l2 / l;
+      if (CGAL::assign(p, result)) {
+        IkFT l = CGAL::approximate_sqrt(sp.data().original_vectors[idx].squared_length());
+
+        IkFT l2 = from_exact(CGAL::approximate_sqrt((p - sp.data().original_rays[idx].point(0)).squared_length()));
+
+        IkFT l3 = (p - sp.data().original_rays[idx].point(0)) * sp.data().original_rays[idx].to_vector();
+        time[i] = l3;
         CGAL_assertion(0 <= time[i]);
-        intersections[i] = from_exact(p);
-        intersections_bary[i] = abs(((from_exact(p) - s) * segment)) / segment_length;
+        intersections[i] = p;
+        intersections_bary[i] = abs(((p - s2) * segment2)) / segment_length2;
         if (!ccw)
           intersections_bary[i] = 1.0 - intersections_bary[i];
       }
@@ -467,12 +488,12 @@ public:
     }
 
     // Calculate pedge vs ivertex collision
-    FT edge_time[2];
+    IkFT edge_time[2];
 
     // Source edge time
     std::size_t adjacent = (source_idx + sp.data().original_vertices.size() - 1) % sp.data().original_vertices.size();
     Vector_2 dir = sp.data().original_vertices[source_idx] - sp.data().original_vertices[adjacent];
-    dir = dir / CGAL::sqrt(dir * dir);
+    dir = dir / CGAL::approximate_sqrt(dir * dir);
 
     // Orthogonal direction matching the direction of the adjacent vertices
     dir = Vector_2(dir.y(), -dir.x());
@@ -491,7 +512,7 @@ public:
     // Target edge time
     adjacent = (target_idx + sp.data().original_vertices.size() - 1) % sp.data().original_vertices.size();
     dir = sp.data().original_vertices[target_idx] - sp.data().original_vertices[adjacent];
-    dir = dir / CGAL::sqrt(dir * dir);
+    dir = dir / CGAL::approximate_sqrt(dir * dir);
 
     // Orthogonal direction matching the direction of the adjacent vertices
     dir = Vector_2(dir.y(), -dir.x());
@@ -509,15 +530,15 @@ public:
 
     // Fill event structure and kinetic intervals.
     if (ccw)
-      kinetic_interval.push_back(std::pair<FT, FT>(0, edge_time[0]));
+      kinetic_interval.push_back(std::pair<IkFT, IkFT>(0, edge_time[0]));
     else
-      kinetic_interval.push_back(std::pair<FT, FT>(0, edge_time[1]));
+      kinetic_interval.push_back(std::pair<IkFT, IkFT>(0, edge_time[1]));
 
     event.time = kinetic_interval.back().second;
     event.intersection_bary = 0;
 
     for (std::size_t i = 0; i < num; i++) {
-      kinetic_interval.push_back(std::pair<FT, FT>(intersections_bary[i], time[i]));
+      kinetic_interval.push_back(std::pair<IkFT, IkFT>(intersections_bary[i], time[i]));
       if (event.time > time[i]) {
         event.time = time[i];
         event.intersection_bary = intersections_bary[i];
@@ -525,23 +546,15 @@ public:
     }
 
     if (ccw)
-      kinetic_interval.push_back(std::pair<FT, FT>(1, edge_time[1]));
+      kinetic_interval.push_back(std::pair<IkFT, IkFT>(1, edge_time[1]));
     else
-      kinetic_interval.push_back(std::pair<FT, FT>(1, edge_time[0]));
+      kinetic_interval.push_back(std::pair<IkFT, IkFT>(1, edge_time[0]));
 
     if (event.time > kinetic_interval.back().second) {
       event.time = kinetic_interval.back().second;
       event.intersection_bary = 1;
     }
 
-/*
-    if (kinetic_interval.size() > 4) {
-      if (kinetic_interval[2].first > kinetic_interval[1].first) {
-        int a;
-        a = 2;
-      }
-    }*/
-
     CGAL_assertion(0 <= event.intersection_bary && event.intersection_bary <= 1);
 
     return event.time;
@@ -549,6 +562,11 @@ public:
 
   template<typename Queue>
   void fill_event_queue(Queue& queue) {
+    // Count faces
+    std::size_t faces = 0;
+    for (std::size_t i = 0; i < m_support_planes.size(); i++)
+      faces += m_support_planes[i].data().mesh.number_of_faces();
+
     for (std::size_t sp_idx = 6; sp_idx < m_support_planes.size(); sp_idx++) {
       std::vector<IEdge> border;
       m_support_planes[sp_idx].get_border(m_intersection_graph, border);
@@ -558,9 +576,10 @@ public:
           continue;
 
         Face_event fe;
-        FT t = calculate_edge_intersection_time(sp_idx, edge, fe);
-        if (t > 0)
+        IkFT t = calculate_edge_intersection_time(sp_idx, edge, fe);
+        if (t > 0) {
           queue.push(fe);
+        }
       }
     }
   }
@@ -746,7 +765,7 @@ public:
 
     remove_equal_points(polygon, 0);
 
-    CGAL_assertion(is_valid_polygon(sp_idx, polygon));
+    is_valid_polygon(sp_idx, polygon);
 
     // Find common planes.
     std::vector<IVertex> vertices;
@@ -989,14 +1008,15 @@ public:
 
   template<typename Pair>
   void sort_points_by_direction(std::vector<Pair>& points) const {
-    FT x = FT(0), y = FT(0);
+    FT x = FT(0), y = FT(0); FT num = 0;
     for (const auto& pair : points) {
       const auto& point = pair.first;
       x += point.x();
       y += point.y();
+      num += 1;
     }
-    x /= static_cast<FT>(points.size());
-    y /= static_cast<FT>(points.size());
+    x /= num;
+    y /= num;
     const Point_2 mid(x, y);
 
     std::sort(points.begin(), points.end(),
@@ -1262,7 +1282,7 @@ public:
     std::size_t line_idx = m_intersection_graph.add_line(line);
     for (std::size_t i = 0; i < vertices.size() - 1; ++i) {
 
-      CGAL_assertion(!is_zero_length_iedge(vertices[i], vertices[i + 1]));
+      //CGAL_assertion(!is_zero_length_iedge(vertices[i], vertices[i + 1]));
       const auto pair = m_intersection_graph.add_edge(
         vertices[i], vertices[i + 1], support_planes_idx);
       const auto iedge = pair.first;
@@ -1381,17 +1401,19 @@ public:
     return support_plane(support_plane_idx).to_2d(segment_3);
   }
 
+/*
   IkSegment_2 to_2d(const std::size_t support_plane_idx, const IkSegment_3& segment_3) const {
     return support_plane(support_plane_idx).to_2d(segment_3);
-  }
+  }*/
 
   Point_2 to_2d(const std::size_t support_plane_idx, const Point_3& point_3) const {
     return support_plane(support_plane_idx).to_2d(point_3);
   }
 
+/*
   IkPoint_2 to_2d(const std::size_t support_plane_idx, const IkPoint_3& point_3) const {
     return support_plane(support_plane_idx).to_2d(point_3);
-  }
+  }*/
 
   Point_2 point_2(const PVertex& pvertex) const {
     return support_plane(pvertex).point_2(pvertex.second);
@@ -1409,9 +1431,10 @@ public:
     return support_plane(support_plane_idx).to_3d(point_2);
   }
 
+/*
   IkPoint_3 to_3d(const std::size_t support_plane_idx, const IkPoint_2& point_2) const {
     return support_plane(support_plane_idx).to_3d(point_2);
-  }
+  }*/
 
   Point_3 point_3(const PVertex& pvertex) const {
     return support_plane(pvertex).point_3(pvertex.second);
@@ -1439,7 +1462,7 @@ public:
     polygon.reserve(points.size());
     for (const auto& pair : points) {
       const auto& p = pair.first;
-      const auto q = to_2d(sp_idx, p);
+      const auto q = m_support_planes[sp_idx].to_2d(p);
       polygon.push_back(std::make_pair(q, true));
     }
     CGAL_assertion(polygon.size() == points.size());
@@ -1450,7 +1473,7 @@ public:
 
     if (!is_valid) {
       for (const auto& pair : polygon) {
-        std::cout << to_3d(sp_idx, pair.first) << std::endl;
+        std::cout << m_support_planes[sp_idx].to_3d(pair.first) << std::endl;
       }
     }
 
diff --git a/Kinetic_space_partition/include/CGAL/KSP_3/Finalizer.h b/Kinetic_space_partition/include/CGAL/KSP_3/Finalizer.h
index f2b4a132db7..44afa6ab175 100644
--- a/Kinetic_space_partition/include/CGAL/KSP_3/Finalizer.h
+++ b/Kinetic_space_partition/include/CGAL/KSP_3/Finalizer.h
@@ -117,6 +117,11 @@ public:
 
     create_volumes();
 
+    if (m_parameters.debug) {
+      for (const auto& v : m_data.volumes())
+        dump_volume(m_data, v.pfaces, "volumes/" + m_data.prefix() + std::to_string(v.index), true, v.index);
+    }
+
     CGAL_assertion(m_data.check_faces());
   }
 
@@ -131,13 +136,15 @@ private:
   void calculate_centroid(Volume_cell& volume) {
     // First find a point in the interior of the volume cell.
     FT x = 0, y = 0, z = 0;
+    FT num = 0;
     for (const PVertex& v : volume.pvertices) {
       Point_3 p = m_data.point_3(v);
       x += p.x();
       y += p.y();
       z += p.z();
+      num += 1;
     }
-    Point_3 inside(x / volume.pvertices.size(), y / volume.pvertices.size(), z / volume.pvertices.size());
+    Point_3 inside(x / num, y / num, z / num);
 
     // Now create a vector of tetrahedrons.
     std::vector<Tetrahedron_3> tets;
@@ -288,7 +295,7 @@ private:
       m_data.incident_faces(m_data.iedge(pedge), neighbor_faces);
 
       if (neighbor_faces.size() == 2) {
-        // If there is only one neighbor, the edge is on the corner of the bbox.
+        // If there is only one neighbor, the edge is on the edge of the bbox.
         // Thus the only neighbor needs to be a bbox face.
         PFace neighbor = (neighbor_faces[0] == pface) ? neighbor_faces[1] : neighbor_faces[0];
         CGAL_assertion(neighbor.first < 6 && pface.first < 6);
@@ -514,21 +521,17 @@ private:
   }
 
   Oriented_side oriented_side(const PFace& a, const PFace& b) const {
-    FT max_dist = 0;
     if (a.first == b.first)
       return CGAL::ON_ORIENTED_BOUNDARY;
     Oriented_side side = CGAL::ON_ORIENTED_BOUNDARY;
-    const Plane_3& p = m_data.support_plane(a.first).plane();
+    const typename Intersection_kernel::Plane_3& p = m_data.support_plane(a.first).exact_plane();
     for (auto v : m_data.pvertices_of_pface(b)) {
-      Point_3 pt = m_data.point_3(v);
-      FT dist = (p.point() - pt) * p.orthogonal_vector();
-      if (CGAL::abs(dist) > max_dist) {
-        side = p.oriented_side(m_data.point_3(v));
-        max_dist = CGAL::abs(dist);
-      }
+      side = p.oriented_side(m_data.point_3(m_data.ivertex(v)));
+      if (side != CGAL::ON_ORIENTED_BOUNDARY)
+        return side;
     }
 
-    return side;
+    return CGAL::ON_ORIENTED_BOUNDARY;
   }
 
   void merge_facets_connected_components() {
diff --git a/Kinetic_space_partition/include/CGAL/KSP_3/Initializer.h b/Kinetic_space_partition/include/CGAL/KSP_3/Initializer.h
index 92b2e92ae33..be34da2cb05 100644
--- a/Kinetic_space_partition/include/CGAL/KSP_3/Initializer.h
+++ b/Kinetic_space_partition/include/CGAL/KSP_3/Initializer.h
@@ -131,6 +131,10 @@ public:
       for (std::size_t sp = 0; sp < m_data.number_of_support_planes(); sp++)
         dump_2d_surface_mesh(m_data, sp, m_data.prefix() + "before-partition-sp" + std::to_string(sp));
     }
+
+    if (m_parameters.verbose) {
+      std::cout << "v: " << m_data.igraph().number_of_vertices() << " f: " << m_data.igraph().number_of_faces() << std::endl;
+    }
   }
 
   void clear() {
@@ -369,7 +373,7 @@ private:
 
   void initial_polygon_iedge_intersections() {
     To_exact to_exact;
-    From_exact to_inexact;
+    From_exact from_exact;
 
     for (std::size_t sp_idx = 0; sp_idx < m_data.number_of_support_planes(); sp_idx++) {
       bool polygons_assigned = false;
@@ -398,11 +402,11 @@ private:
         typename Intersection_kernel::Point_2 a(sp.to_2d(m_data.point_3(m_data.source(pair.second[0]))));
         typename Intersection_kernel::Point_2 b(sp.to_2d(m_data.point_3(m_data.target(pair.second[0]))));
         typename Intersection_kernel::Line_2 exact_line(a, b);
-        Line_2 l = to_inexact(exact_line);
+        Line_2 l = from_exact(exact_line);
 
         typename Intersection_kernel::Vector_2 ldir = exact_line.to_vector();
         ldir = (typename Intersection_kernel::FT(1.0) / CGAL::approximate_sqrt(ldir * ldir)) * ldir;
-        Vector_2 dir = to_inexact(ldir);
+        Vector_2 dir = from_exact(ldir);
 
         std::vector<typename Intersection_kernel::Segment_2> crossing_polygon_segments;
         std::vector<IEdge> crossing_iedges;
@@ -433,20 +437,20 @@ private:
               if (eproj < emin) {
                 eminp = intersection;
                 emin = eproj;
-                minp = to_inexact(intersection);
+                minp = from_exact(intersection);
                 //min = proj;
                 typename Intersection_kernel::FT p = dir * edge_dir;
                 assert(p != 0);
-                min_speed = CGAL::sqrt(edge_dir * edge_dir) / to_inexact(p);
+                min_speed = CGAL::approximate_sqrt(edge_dir * edge_dir) / from_exact(p);
               }
               if (emax < eproj) {
                 emaxp = intersection;
                 emax = eproj;
-                maxp = to_inexact(intersection);
+                maxp = from_exact(intersection);
                 //max = proj;
                 typename Intersection_kernel::FT p = dir * edge_dir;
                 assert(p != 0);
-                max_speed = CGAL::sqrt(edge_dir * edge_dir) / to_inexact(p);
+                max_speed = CGAL::approximate_sqrt(edge_dir * edge_dir) / from_exact(p);
               }
             }
             else std::cout << "crossing segment does not intersect line" << std::endl;
@@ -497,9 +501,9 @@ private:
                 crossing_iedges.push_back(e);
                 if (emin > s) {
                   typename Intersection_kernel::FT bary_edge_exact = (emin - s) / (t - s);
-                  FT bary_edge = to_inexact((emin - s) / (t - s));
+                  FT bary_edge = from_exact((emin - s) / (t - s));
                   CGAL_assertion(bary_edge_exact >= 0);
-                  FT time = CGAL::abs(to_inexact(s - emin) / min_speed);
+                  FT time = CGAL::abs(from_exact(s - emin) / min_speed);
                   kinetic_interval.push_back(std::pair<FT, FT>(0, time)); // border barycentric coordinate
                   kinetic_interval.push_back(std::pair<FT, FT>(bary_edge, 0));
                 }
@@ -509,9 +513,9 @@ private:
 
                 if (t > emax) {
                   typename Intersection_kernel::FT bary_edge_exact = (emax - s) / (t - s);
-                  FT bary_edge = to_inexact((emax - s) / (t - s));
+                  FT bary_edge = from_exact((emax - s) / (t - s));
                   CGAL_assertion(0 <= bary_edge_exact && bary_edge_exact <= 1);
-                  FT time = CGAL::abs(to_inexact(emax - t) / max_speed);
+                  FT time = CGAL::abs(from_exact(emax - t) / max_speed);
                   kinetic_interval.push_back(std::pair<FT, FT>(bary_edge, 0));
                   kinetic_interval.push_back(std::pair<FT, FT>(1, time)); // border barycentric coordinate
                 }
@@ -541,9 +545,9 @@ private:
               crossing_iedges.push_back(e);
               if (s > emax) {
                 typename Intersection_kernel::FT bary_edge_exact = (s - emax) / (s - t);
-                FT bary_edge = to_inexact((s - emax) / (s - t));
+                FT bary_edge = from_exact((s - emax) / (s - t));
                 CGAL_assertion(0 <= bary_edge_exact && bary_edge_exact <= 1);
-                FT time = CGAL::abs(to_inexact(emax - s) / max_speed);
+                FT time = CGAL::abs(from_exact(emax - s) / max_speed);
                 kinetic_interval.push_back(std::pair<FT, FT>(0, time)); // border barycentric coordinate
                 kinetic_interval.push_back(std::pair<FT, FT>(bary_edge, 0));
               }
@@ -552,9 +556,9 @@ private:
 
               if (emin > t) {
                 typename Intersection_kernel::FT bary_edge_exact = (s - emin) / (s - t);
-                FT bary_edge = to_inexact(bary_edge_exact);
+                FT bary_edge = from_exact(bary_edge_exact);
                 CGAL_assertion(0 <= bary_edge_exact && bary_edge_exact <= 1);
-                FT time = CGAL::abs(to_inexact(t - emin) / min_speed);
+                FT time = CGAL::abs(from_exact(t - emin) / min_speed);
                 kinetic_interval.push_back(std::pair<FT, FT>(bary_edge, 0));
                 kinetic_interval.push_back(std::pair<FT, FT>(1, time)); // border barycentric coordinate
               }
@@ -832,11 +836,11 @@ private:
           typename Intersection_kernel::Point_2 point;
           typename Intersection_kernel::Segment_3 seg_a(m_data.point_3(it_a->second.first), m_data.point_3(it_a->second.second));
           typename Intersection_kernel::Segment_3 seg_b(m_data.point_3(it_b->second.first), m_data.point_3(it_b->second.second));
-          if (!intersection(m_data.to_2d(common_plane_idx, seg_a), m_data.to_2d(common_plane_idx, seg_b), point))
+          if (!intersection(m_data.support_plane(common_plane_idx).to_2d(seg_a), m_data.support_plane(common_plane_idx).to_2d(seg_b), point))
             continue;
 
           crossed_vertices.push_back(
-            m_data.add_ivertex(m_data.to_3d(common_plane_idx, point), union_set));
+            m_data.add_ivertex(m_data.support_plane(common_plane_idx).to_3d(point), union_set));
         }
       }
       crossed_vertices.push_back(it_a->second.second);
diff --git a/Kinetic_space_partition/include/CGAL/KSP_3/Intersection_graph.h b/Kinetic_space_partition/include/CGAL/KSP_3/Intersection_graph.h
index 207ea0d2c0f..e86d696bc0b 100644
--- a/Kinetic_space_partition/include/CGAL/KSP_3/Intersection_graph.h
+++ b/Kinetic_space_partition/include/CGAL/KSP_3/Intersection_graph.h
@@ -39,21 +39,20 @@ public:
   using Kernel = GeomTraits;
   using Intersection_kernel = IntersectionKernel;
 
+  using IkFT = typename Intersection_kernel::FT;
   using Point_2 = typename Intersection_kernel::Point_2;
   using Point_3 = typename Intersection_kernel::Point_3;
   using Segment_3 = typename Intersection_kernel::Segment_3;
   using Line_3 = typename Intersection_kernel::Line_3;
   using Polygon_2 = typename CGAL::Polygon_2<Intersection_kernel>;
 
-  using Inexact_FT = typename Kernel::FT;
-
   struct Vertex_property {
     Point_3 point;
     Vertex_property() {}
     Vertex_property(const Point_3& point) : point(point) {}
   };
 
-  using Kinetic_interval = std::vector<std::pair<Inexact_FT, Inexact_FT> >;
+  using Kinetic_interval = std::vector<std::pair<IkFT, IkFT> >;
 
   struct Edge_property {
     std::size_t line;
diff --git a/Kinetic_space_partition/include/CGAL/KSP_3/FacePropagation.h b/Kinetic_space_partition/include/CGAL/KSP_3/Propagation.h
similarity index 92%
rename from Kinetic_space_partition/include/CGAL/KSP_3/FacePropagation.h
rename to Kinetic_space_partition/include/CGAL/KSP_3/Propagation.h
index c6a2816b3da..26e9e228a8a 100644
--- a/Kinetic_space_partition/include/CGAL/KSP_3/FacePropagation.h
+++ b/Kinetic_space_partition/include/CGAL/KSP_3/Propagation.h
@@ -30,7 +30,7 @@ namespace internal {
 #else
 
 template<typename GeomTraits, typename IntersectionKernel>
-class FacePropagation {
+class Propagation {
 
 public:
   using Kernel = GeomTraits;
@@ -38,6 +38,7 @@ public:
 
 private:
   using FT = typename Kernel::FT;
+  using IkFT = typename Intersection_kernel::FT;
   using Point_2 = typename Kernel::Point_2;
   using Vector_2 = typename Kernel::Vector_2;
   using Segment_2 = typename Kernel::Segment_2;
@@ -70,7 +71,7 @@ private:
   };
 
 public:
-  FacePropagation(Data_structure& data, const Parameters& parameters) :
+  Propagation(Data_structure& data, const Parameters& parameters) :
     m_data(data), m_parameters(parameters),
     m_min_time(-FT(1)), m_max_time(-FT(1))
   { }
@@ -182,8 +183,8 @@ private:
 
           // Within an interval
           if (ki->second[i].first > event.intersection_bary && ki->second[i - 1].first < event.intersection_bary) {
-            FT interval_pos = (event.intersection_bary - ki->second[i - 1].first) / (ki->second[i].first - ki->second[i - 1].first);
-            FT interval_time = interval_pos * (ki->second[i].second - ki->second[i - 1].second) + ki->second[i - 1].second;
+            IkFT interval_pos = (event.intersection_bary - ki->second[i - 1].first) / (ki->second[i].first - ki->second[i - 1].first);
+            IkFT interval_time = interval_pos * (ki->second[i].second - ki->second[i - 1].second) + ki->second[i - 1].second;
 
             if (event.time > interval_time) {
               crossing++;
@@ -215,7 +216,7 @@ private:
 
     for (IEdge edge : border) {
       Face_event fe;
-      FT t = m_data.calculate_edge_intersection_time(event.support_plane, edge, fe);
+      IkFT t = m_data.calculate_edge_intersection_time(event.support_plane, edge, fe);
       if (t > 0)
         m_face_queue.push(fe);
     }
diff --git a/Kinetic_space_partition/include/CGAL/KSP_3/Support_plane.h b/Kinetic_space_partition/include/CGAL/KSP_3/Support_plane.h
index dba5675bbe1..d719f56307a 100644
--- a/Kinetic_space_partition/include/CGAL/KSP_3/Support_plane.h
+++ b/Kinetic_space_partition/include/CGAL/KSP_3/Support_plane.h
@@ -82,8 +82,8 @@ public:
     Face_event() {}
     Face_event(std::size_t sp_idx, FT time, IEdge edge, IFace face) : support_plane(sp_idx), time(time), crossed_edge(edge), face(face) {}
     std::size_t support_plane;
-    FT time;
-    FT intersection_bary;
+    typename Intersection_kernel::FT time;
+    typename Intersection_kernel::FT intersection_bary;
     IEdge crossed_edge;
     IFace face; // The face that does not yet belong to the region.
   };
@@ -132,6 +132,7 @@ public:
   };
 
 private:
+  static constexpr bool identical_kernel = !std::is_same_v<Kernel, Intersection_kernel>;
 
   std::shared_ptr<Data> m_data;
 
@@ -172,55 +173,9 @@ public:
     add_property_maps();
   }
 
-  template<typename PointRange>
-  Support_plane(const PointRange& polygon, const bool is_bbox) :
-    m_data(std::make_shared<Data>()) {
-    To_exact to_exact;
-
-    std::vector<Point_3> points;
-    points.reserve(polygon.size());
-    for (const auto& point : polygon) {
-      points.push_back(Point_3(
-        static_cast<FT>(point.x()),
-        static_cast<FT>(point.y()),
-        static_cast<FT>(point.z())));
-    }
-    const std::size_t n = points.size();
-    CGAL_assertion(n == polygon.size());
-
-    Vector_3 normal = CGAL::NULL_VECTOR;
-    for (std::size_t i = 0; i < n; ++i) {
-      const std::size_t ip = (i + 1) % n;
-      const auto& pa = points[i];
-      const auto& pb = points[ip];
-      const FT x = normal.x() + (pa.y() - pb.y()) * (pa.z() + pb.z());
-      const FT y = normal.y() + (pa.z() - pb.z()) * (pa.x() + pb.x());
-      const FT z = normal.z() + (pa.x() - pb.x()) * (pa.y() + pb.y());
-      normal = Vector_3(x, y, z);
-    }
-    CGAL_assertion_msg(normal != CGAL::NULL_VECTOR, "ERROR: BBOX IS FLAT!");
-    CGAL_assertion(n != 0);
-
-    m_data->k = 0;
-    m_data->plane = Plane_3(points[0], KSP::internal::normalize(normal));
-    m_data->exact_plane = to_exact(m_data->plane);
-    m_data->is_bbox = is_bbox;
-    m_data->distance_tolerance = 0;
-    m_data->angle_tolerance = 0;
-    m_data->actual_input_polygon = -1;
-
-    std::vector<Triangle_2> tris(points.size() - 2);
-    for (std::size_t i = 2; i < points.size(); i++) {
-      tris[i - 2] = Triangle_2(to_2d(points[0]), to_2d(points[i - 1]), to_2d(points[i]));
-    }
-
-    m_data->centroid = CGAL::centroid(tris.begin(), tris.end(), CGAL::Dimension_tag<2>());
-
-    add_property_maps();
-  }
-
   Support_plane(const std::vector<typename Intersection_kernel::Point_3>& polygon, const bool is_bbox) :
     m_data(std::make_shared<Data>()) {
+
     From_exact from_exact;
 
     std::vector<Point_3> points;
@@ -413,15 +368,17 @@ public:
 
     std::vector<std::pair<std::size_t, Direction_2> > dir_vec;
 
+    FT num = 0;
     for (const auto& pair : points) {
       const auto& point = pair.first;
       directions.push_back(Vector_2(m_data->centroid, point));
       const FT length = static_cast<FT>(
-        CGAL::sqrt(CGAL::to_double(CGAL::abs(directions.back() * directions.back()))));
+        CGAL::approximate_sqrt(CGAL::abs(directions.back() * directions.back())));
       sum_length += length;
+      num += 1;
     }
     CGAL_assertion(directions.size() == n);
-    sum_length /= static_cast<FT>(n);
+    sum_length /= num;
 
     dir_vec.reserve(n);
     for (std::size_t i = 0; i < n; i++)
@@ -439,7 +396,7 @@ public:
       m_data->original_vertices[i] = point;
       m_data->original_vectors[i] = directions[dir_vec[i].first] / sum_length;
       m_data->original_directions[i] = Direction_2(directions[dir_vec[i].first]);
-      m_data->original_rays[i] = typename Intersection_kernel::Ray_2(to_exact(point), to_exact(m_data->original_directions[i]));
+      m_data->original_rays[i] = typename Intersection_kernel::Ray_2(to_exact(point), to_exact(m_data->original_vectors[i]));
       m_data->v_original_map[vi] = true;
       vertices.push_back(vi);
     }
@@ -667,7 +624,7 @@ public:
   const Vector_2 original_edge_direction(std::size_t v1, std::size_t v2) const {
     const Vector_2 edge = m_data->original_vertices[v1] - m_data->original_vertices[v2];
     Vector_2 orth = Vector_2(-edge.y(), edge.x());
-    orth = (1.0 / (CGAL::sqrt(orth * orth))) * orth;
+    orth = (1.0 / (CGAL::approximate_sqrt(orth * orth))) * orth;
     FT s1 = orth * m_data->original_vectors[v1];
     FT s2 = orth * m_data->original_vectors[v2];
 
@@ -678,8 +635,7 @@ public:
   }
 
   const FT speed(const Vertex_index& vi) const {
-    return static_cast<FT>(CGAL::sqrt(
-      CGAL::to_double(CGAL::abs(m_data->direction[vi].squared_length()))));
+    return static_cast<FT>(CGAL::approximate_sqrt((CGAL::abs(m_data->direction[vi].squared_length()))));
   }
 
   const std::vector<std::size_t>& input(const Face_index& fi) const { return m_data->input_map[fi]; }
@@ -712,6 +668,7 @@ public:
       m_data->plane.to_2d(Point_3(0, 0, 0) + vec));
   }
 
+  template<typename = typename std::enable_if<identical_kernel>::type >
   const typename Intersection_kernel::Point_2 to_2d(const typename Intersection_kernel::Point_3& point) const {
     return m_data->exact_plane.to_2d(point);
   }
@@ -722,6 +679,7 @@ public:
       m_data->plane.to_2d(line.point() + line.to_vector()));
   }
 
+  template<typename = typename std::enable_if<identical_kernel>::type >
   const typename Intersection_kernel::Line_2 to_2d(const typename Intersection_kernel::Line_3& line) const {
     return typename Intersection_kernel::Line_2(
       m_data->exact_plane.to_2d(line.point()),
@@ -734,6 +692,7 @@ public:
       m_data->plane.to_2d(segment.target()));
   }
 
+  template<typename = typename std::enable_if<identical_kernel>::type >
   const typename Intersection_kernel::Segment_2 to_2d(const typename Intersection_kernel::Segment_3& segment) const {
     return typename Intersection_kernel::Segment_2(
       m_data->exact_plane.to_2d(segment.source()),
@@ -750,6 +709,7 @@ public:
     return m_data->plane.to_3d(point);
   }
 
+  template<typename = typename std::enable_if<identical_kernel>::type >
   const typename Intersection_kernel::Point_3 to_3d(const typename Intersection_kernel::Point_2& point) const {
     return m_data->exact_plane.to_3d(point);
   }
@@ -769,18 +729,6 @@ public:
   const Vertex_index add_vertex(const Point_2& point) {
     return m_data->mesh.add_vertex(point);
   }
-/*
-
-  const Vertex_index duplicate_vertex(const Vertex_index& v) {
-    // TODO: We cannot take it by reference because it fails for EPECK
-    // when called from  front_and_back_34() in Data_structure.
-    const auto pp = m_data->mesh.point(v);
-    const auto vi = m_data->mesh.add_vertex(pp);
-    m_data->direction[vi] = m_data->direction[v];
-    m_data->v_ivertex_map[vi] = m_data->v_ivertex_map[v];
-    m_data->v_iedge_map[vi] = m_data->v_iedge_map[v];
-    return vi;
-  }*/
 
   void remove_vertex(const Vertex_index& vi) {
     m_data->mesh.remove_vertex(vi);
@@ -809,38 +757,9 @@ bool operator==(const Support_plane<GeomTraits, IntersectionKernel>& a, const Su
     return false;
   }
 
-  using FT = typename GeomTraits::FT;
-  const auto& planea = a.plane();
-  const auto& planeb = b.plane();
-
-  const auto va = planea.orthogonal_vector();
-  const auto vb = planeb.orthogonal_vector();
-
-  // Are the planes parallel?
-
-  FT aval = approximate_angle(va, vb);
-  CGAL_assertion(aval >= FT(0) && aval <= FT(180));
-  if (aval >= FT(90))
-    aval = FT(180) - aval;
-
-  if (aval >= a.angle_tolerance()) {
-    return false;
-  }
-
-  const auto pa1 = a.to_3d(a.centroid());
-  const auto pb1 = planeb.projection(pa1);
-  const auto pb2 = b.to_3d(b.centroid());
-  const auto pa2 = planea.projection(pb2);
-
-  const FT bval1 = KSP::internal::distance(pa1, pb1);
-  const FT bval2 = KSP::internal::distance(pa2, pb2);
-  const FT bval = (CGAL::max)(bval1, bval2);
-  CGAL_assertion(bval >= FT(0));
-
-  if (bval >= a.distance_tolerance())
-    return false;
-
-  return true;
+  if (a.exact_plane() == b.exact_plane())
+    return true;
+  else return false;
 }
 
 #endif //DOXYGEN_RUNNING
diff --git a/Kinetic_space_partition/include/CGAL/Kinetic_space_partition_3.h b/Kinetic_space_partition/include/CGAL/Kinetic_space_partition_3.h
index da79f6f51ad..4de442bd399 100644
--- a/Kinetic_space_partition/include/CGAL/Kinetic_space_partition_3.h
+++ b/Kinetic_space_partition/include/CGAL/Kinetic_space_partition_3.h
@@ -18,7 +18,6 @@
 // Boost includes.
 #include <CGAL/boost/graph/named_params_helper.h>
 #include <CGAL/Named_function_parameters.h>
-//#include <boost/filesystem.hpp>
 
 #include <algorithm>
 #include <numeric>
@@ -44,7 +43,7 @@
 
 #include <CGAL/KSP_3/Data_structure.h>
 #include <CGAL/KSP_3/Initializer.h>
-#include <CGAL/KSP_3/FacePropagation.h>
+#include <CGAL/KSP_3/Propagation.h>
 #include <CGAL/KSP_3/Finalizer.h>
 
 #include <CGAL/Octree.h>
@@ -136,7 +135,7 @@ private:
   using To_exact = typename CGAL::Cartesian_converter<Kernel, Intersection_kernel>;
 
   using Initializer = KSP_3::internal::Initializer<Kernel, Intersection_kernel>;
-  using Propagation = KSP_3::internal::FacePropagation<Kernel, Intersection_kernel>;
+  using Propagation = KSP_3::internal::Propagation<Kernel, Intersection_kernel>;
   using Finalizer   = KSP_3::internal::Finalizer<Kernel, Intersection_kernel>;
 
   using Polygon_mesh = CGAL::Surface_mesh<Point_3>;
@@ -482,7 +481,6 @@ public:
     typename NamedParameters = parameters::Default_named_parameters>
   void initialize(
     const NamedParameters& np = CGAL::parameters::default_values()) {
-
     Timer timer;
     m_parameters.bbox_dilation_ratio = parameters::choose_parameter(
       parameters::get_parameter(np, internal_np::bbox_dilation_ratio), FT(11) / FT(10));
@@ -529,7 +527,7 @@ public:
 
     if (m_parameters.debug) {
       for (std::size_t i = 0; i < m_input_polygons.size(); i++)
-        KSP_3::internal::dump_polygon(m_input_polygons[i], std::to_string(i) + "-input_polygon");
+        KSP_3::internal::dump_polygon<GeomTraits>(m_input_polygons[i], std::to_string(i) + "-input_polygon");
     }
 
     split_octree();
@@ -708,38 +706,13 @@ public:
 
     timer.stop();
 
- /*   if (m_parameters.debug) {
-
-      if (boost::filesystem::is_directory("volumes/"))
-        for (boost::filesystem::directory_iterator end_dir_it, it("volumes/"); it != end_dir_it; ++it)
-          boost::filesystem::remove_all(it->path());
-
-      KSP_3::dump_volumes_ksp(*this, "volumes/");
-      for (std::size_t i = 1; i < m_volumes.size(); i++)
-        if (m_volumes[i].first != m_volumes[i - 1].first)
-          std::cout << i << " " << m_volumes[i - 1].first << std::endl;
-      std::cout << m_volumes.size() << " " << m_volumes.back().first << std::endl;
-    }*/
-
     timer.reset();
     timer.start();
     make_conformal(0);
     conformal_time = timer.time();
 
-/*    if (m_parameters.debug) {
-
-      if (boost::filesystem::is_directory("volumes_after/"))
-        for (boost::filesystem::directory_iterator end_dir_it, it("volumes_after/"); it != end_dir_it; ++it)
-          boost::filesystem::remove_all(it->path());
-      KSP_3::dump_volumes_ksp(*this, "volumes_after/");
-      for (std::size_t i = 1; i < m_volumes.size(); i++)
-        if (m_volumes[i].first != m_volumes[i - 1].first)
-          std::cout << i << " " << m_volumes[i - 1].first << std::endl;
-      std::cout << m_volumes.size() << " " << m_volumes.back().first << std::endl;
-    }*/
-
-    if (m_parameters.verbose)
-      check_tjunctions();
+//     if (m_parameters.verbose)
+//       check_tjunctions();
 
     // Clear unused data structures
     for (std::size_t i = 0; i < m_partitions.size(); i++) {
@@ -749,6 +722,8 @@ public:
       m_partition_nodes[i].m_data->face_to_volumes().clear();
     }
 
+    std::cout << "ksp v: " << m_partition_nodes[0].m_data->vertices().size() << " f: " << m_partition_nodes[0].face2vertices.size() << " vol: " << m_volumes.size() << std::endl;
+
     return;
   }
 #endif
@@ -802,7 +777,12 @@ public:
     std::vector<typename Intersection_kernel::Point_3> vtx;
     std::vector<Index> vtx_index;
 
-    From_exact to_inexact;
+    using LCC_kernel = typename CGAL::Kernel_traits<typename LCC::Point>::Kernel;
+
+    using To_lcc = CGAL::Cartesian_converter<Intersection_kernel, LCC_kernel>;
+
+    To_lcc to_lcc;
+
     To_exact to_exact;
 
     std::vector<Index> faces_of_volume, vtx_of_face;
@@ -832,7 +812,7 @@ public:
 
     CGAL::Linear_cell_complex_incremental_builder_3<LCC> ib(lcc);
     for (std::size_t i = 0; i < vtx.size(); i++)
-      ib.add_vertex(vtx[i]);
+      ib.add_vertex(to_lcc(vtx[i]));
 
     std::size_t num_faces = 0;
     //std::size_t num_vols = 0;
@@ -905,12 +885,12 @@ public:
           std::size_t nn = (n + 1) % vtx_of_face.size();
           norm = CGAL::cross_product(vtx[mapped_vertices[vtx_of_face[n]]] - vtx[mapped_vertices[vtx_of_face[i]]], vtx[mapped_vertices[vtx_of_face[nn]]] - vtx[mapped_vertices[vtx_of_face[n]]]);
           i++;
-        } while (to_inexact(norm.squared_length()) == 0 && i < vtx_of_face.size());
+        } while (norm.squared_length() == 0 && i < vtx_of_face.size());
 
-        FT len = sqrt(to_inexact(norm.squared_length()));
+        typename Intersection_kernel::FT len = CGAL::approximate_sqrt(norm.squared_length());
         if (len != 0)
           len = 1.0 / len;
-        norm = norm * to_exact(len);
+        norm = norm * len;
 
         bool outwards_oriented = (vtx[mapped_vertices[vtx_of_face[0]]] - centroid) * norm < 0;
         //outward[std::make_pair(v, j)] = outwards_oriented;
@@ -1856,9 +1836,9 @@ private:
 
     Vector_2 axis1 = bbox[0] - bbox[1];
     Vector_2 axis2 = bbox[1] - bbox[2];
-    FT la = CGAL::sqrt(axis1.squared_length());
+    FT la = CGAL::approximate_sqrt(axis1.squared_length());
     axis1 = axis1 * (1.0 / la);
-    FT lb = CGAL::sqrt(axis2.squared_length());
+    FT lb = CGAL::approximate_sqrt(axis2.squared_length());
     axis2 = axis2 * (1.0 / lb);
 
     if (CGAL::abs(axis1.x()) < CGAL::abs(axis2.x())) {
@@ -2294,7 +2274,6 @@ private:
     From_exact from_exact;
 
     if (m_parameters.reorient_bbox) {
-
       m_transform = to_exact(get_obb2abb(m_input_polygons));
 
       for (const auto& p : m_input_polygons) {
@@ -2407,7 +2386,7 @@ private:
           vout << std::endl;
           vout.close();
 
-          KSP_3::internal::dump_polygons(m_partition_nodes[idx].clipped_polygons, std::to_string(idx) + "-polys.ply");
+          //KSP_3::internal::dump_polygons(m_partition_nodes[idx].clipped_polygons, std::to_string(idx) + "-polys.ply");
         }
         idx++;
       }
diff --git a/Kinetic_space_partition/test/Kinetic_space_partition/kinetic_3d_test_all.cpp b/Kinetic_space_partition/test/Kinetic_space_partition/kinetic_3d_test_all.cpp
index 9263cf127b1..77531de7cc0 100644
--- a/Kinetic_space_partition/test/Kinetic_space_partition/kinetic_3d_test_all.cpp
+++ b/Kinetic_space_partition/test/Kinetic_space_partition/kinetic_3d_test_all.cpp
@@ -1,25 +1,25 @@
 #include <CGAL/Simple_cartesian.h>
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/intersections.h>
 #include <CGAL/Kinetic_space_partition_3.h>
-#include <CGAL/Real_timer.h>
 #include <CGAL/IO/OFF.h>
 #include <CGAL/IO/PLY.h>
 
-using SCF   = CGAL::Simple_cartesian<float>;
 using SCD   = CGAL::Simple_cartesian<double>;
 using EPICK = CGAL::Exact_predicates_inexact_constructions_kernel;
 using EPECK = CGAL::Exact_predicates_exact_constructions_kernel;
-using Timer = CGAL::Real_timer;
 
-template<typename Kernel>
+std::size_t different = 0;
+
+template<typename Kernel, typename IntersectionKernel>
 bool run_test(
   const std::string input_filename,
   const std::vector<unsigned int>& ks,
   const std::vector<std::vector<unsigned int> >& results) {
 
   using Point_3   = typename Kernel::Point_3;
-  using KSP = CGAL::Kinetic_space_partition_3<Kernel>;
+  using KSP = CGAL::Kinetic_space_partition_3<Kernel, IntersectionKernel>;
 
   std::string filename = input_filename;
   std::ifstream input_file_off(filename);
@@ -47,7 +47,7 @@ bool run_test(
     //std::cout << std::endl << "--INPUT K: " << k << std::endl;
     ksp.partition(ks[i]);
 
-    CGAL::Linear_cell_complex_for_combinatorial_map<3, 3, CGAL::Linear_cell_complex_traits<3, CGAL::Exact_predicates_exact_constructions_kernel>, typename KSP::Linear_cell_complex_min_items> lcc;
+    CGAL::Linear_cell_complex_for_combinatorial_map<3, 3, CGAL::Linear_cell_complex_traits<3, Kernel>, typename KSP::Linear_cell_complex_min_items> lcc;
     ksp.get_linear_cell_complex(lcc);
 
     std::vector<unsigned int> cells = { 0, 2, 3 }, count;
@@ -56,243 +56,237 @@ bool run_test(
     std::cout << ksp.number_of_volumes() << std::endl;
 
     if (results[i][0] != count[0] || results[i][1] != count[2] || results[i][2] != count[3]) {
-      std::cout << "TEST FAILED: Partitioning has not expected number of vertices, faces or volumes for k = " << ks[i] << std::endl;
+      std::cout << "TEST differs: Partitioning has not expected number of vertices, faces or volumes for k = " << ks[i] << std::endl;
 
       std::cout << "Expectation:" << std::endl;
       std::cout << "v: " << results[i][0] << " f : " << results[i][1] << " v : " << results[i][2] << std::endl;
       std::cout << "Result k = " << " vertices : " << count[0] << " faces : " << count[2] << " volumes : " << count[3] << std::endl;
-      //assert(false);
+      std::cout << input_filename << std::endl;
+
+      different++;
     }
+    else std::cout << "TEST PASSED k = " << ks[i] << " " << input_filename << std::endl;
   }
 
   return true;
 }
 
-template<typename Kernel>
+template<typename Kernel, typename IntersectionKernel>
 void run_all_tests() {
+  different = 0;
   std::cout.precision(10);
   std::vector< std::vector<double> > all_times;
 
-  // All results are precomputed for k = 1!
-  std::vector<std::vector<unsigned int> > results(3); //
+  std::vector<std::vector<unsigned int> > results(3);
 
-  //run_test<Kernel>("20-inserted-polygons.ply", { 3 }, results);
-
-  results[0] = { 58, 89, 20 };
-  results[1] = { 63, 102, 24 };
-  results[2] = { 63, 106, 26 };
-  run_test<Kernel>("data/stress-test-5/test-2-rnd-polygons-20-4.off", { 1, 2, 3 }, results);
-
-  results[0] = { 206, 385, 99 };
-  results[1] = { 232, 449, 118 };
-  results[2] = { 265, 540, 147 };
-  run_test<Kernel>("data/real-data-test/test-15-polygons.off", { 1, 2, 3 }, results);
-
-  results[0] = { 39, 49, 10 };
+  results[0] = { 40, 52, 11 };
   results[1] = { 48, 70, 16 };
   results[2] = { 54, 84, 20 };
-  run_test<Kernel>("data/edge-case-test/test-same-time.off", { 1, 2, 3 }, results);
+  run_test<Kernel, IntersectionKernel>("data/edge-case-test/test-same-time.off", { 1, 2, 3 }, results);
 
   // Edge tests.
   results[0] = { 18, 20, 4 };
-  run_test<Kernel>("data/edge-case-test/test-2-polygons.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/edge-case-test/test-2-polygons.off", { 1 }, results);
 
   results[0] = { 22, 25, 5 };
-  run_test<Kernel>("data/edge-case-test/test-4-polygons.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/edge-case-test/test-4-polygons.off", { 1 }, results);
 
   results[0] = { 22, 25, 5 };
-  run_test<Kernel>("data/edge-case-test/test-5-polygons.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/edge-case-test/test-5-polygons.off", { 1 }, results);
 
-  results[0] = { 40, 52, 11 };
-  results[1] = { 51, 77, 18 };
-  run_test<Kernel>("data/edge-case-test/test-local-global-1.off", { 1, 2 }, results);
-
-  results[0] = { 40, 52, 11 };
+  results[0] = { 39, 49, 10 };
   results[1] = { 49, 73, 17 };
+  run_test<Kernel, IntersectionKernel>("data/edge-case-test/test-local-global-1.off", { 1, 2 }, results);
+
+  results[0] = { 39, 49, 10 };
+  results[1] = { 51, 77, 18 };
   results[2] = { 54, 84, 20 };
-  run_test<Kernel>("data/edge-case-test/test-local-global-2.off", { 1, 2, 3 }, results);
+  run_test<Kernel, IntersectionKernel>("data/edge-case-test/test-local-global-2.off", { 1, 2, 3 }, results);
 
   // Stress tests 0.
   results[0] = { 14, 13, 2 };
-  run_test<Kernel>("data/stress-test-0/test-1-polygon-a.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-0/test-1-polygon-a.off", { 1 }, results);
   results[0] = { 14, 13, 2 };
-  run_test<Kernel>("data/stress-test-0/test-1-polygon-b.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-0/test-1-polygon-b.off", { 1 }, results);
   results[0] = { 14, 13, 2 };
-  run_test<Kernel>("data/stress-test-0/test-1-polygon-c.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-0/test-1-polygon-c.off", { 1 }, results);
   results[0] = { 14, 13, 2 };
-  run_test<Kernel>("data/stress-test-0/test-1-polygon-d.off", { 1 }, results);
-  results[0] = { 20, 22, 4 };
-  run_test<Kernel>("data/stress-test-0/test-2-polygons-ab.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-0/test-1-polygon-d.off", { 1 }, results);
+  results[0] = { 18, 18, 3 };
+  run_test<Kernel, IntersectionKernel>("data/stress-test-0/test-2-polygons-ab.off", { 1 }, results);
   results[0] = { 19, 19, 3 };
   results[1] = { 20, 22, 4 };
-  run_test<Kernel>("data/stress-test-0/test-2-polygons-ac.off", { 1, 2 }, results);
-  results[0] = { 20, 22, 4 };
-  run_test<Kernel>("data/stress-test-0/test-2-polygons-ad.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-0/test-2-polygons-ac.off", { 1, 2 }, results);
+  results[0] = { 19, 19, 3 };
+  run_test<Kernel, IntersectionKernel>("data/stress-test-0/test-2-polygons-ad.off", { 1 }, results);
   results[0] = { 18, 18, 3 };
-  run_test<Kernel>("data/stress-test-0/test-2-polygons-bc.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-0/test-2-polygons-bc.off", { 1 }, results);
   results[0] = { 18, 18, 3 };
   results[1] = { 19, 21, 4 };
-  run_test<Kernel>("data/stress-test-0/test-2-polygons-bd.off", { 1, 2 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-0/test-2-polygons-bd.off", { 1, 2 }, results);
   results[0] = { 19, 21, 4 };
-  run_test<Kernel>("data/stress-test-0/test-2-polygons-cd.off", { 1 }, results);
-  results[0] = { 27, 32, 6 };
-  run_test<Kernel>("data/stress-test-0/test-3-polygons-abc.off", { 1 }, results);
-  results[0] = { 28, 33, 6 };
+  run_test<Kernel, IntersectionKernel>("data/stress-test-0/test-2-polygons-cd.off", { 1 }, results);
+  results[0] = { 26, 29, 5 };
+  run_test<Kernel, IntersectionKernel>("data/stress-test-0/test-3-polygons-abc.off", { 1 }, results);
+  results[0] = { 28, 31, 5 };
   results[1] = { 30, 39, 8 };
-  run_test<Kernel>("data/stress-test-0/test-3-polygons-abd.off", { 1, 2 }, results);
-  results[0] = { 27, 32, 6 };
-  results[1] = { 28, 35, 7 };
-  run_test<Kernel>("data/stress-test-0/test-3-polygons-acd.off", { 1, 2 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-0/test-3-polygons-abd.off", { 1, 2 }, results);
+  results[0] = { 25, 28, 5 };
+  results[1] = { 27, 32, 6 };
+  run_test<Kernel, IntersectionKernel>("data/stress-test-0/test-3-polygons-acd.off", { 1, 2 }, results);
   results[0] = { 25, 28, 5 };
   results[1] = { 26, 31, 6 };
-  run_test<Kernel>("data/stress-test-0/test-3-polygons-bcd.off", { 1, 2 }, results);
-  results[0] = { 36, 46, 9 };
+  run_test<Kernel, IntersectionKernel>("data/stress-test-0/test-3-polygons-bcd.off", { 1, 2 }, results);
+  results[0] = { 34, 42, 8 };
   results[1] = { 38, 52, 11 };
-  run_test<Kernel>("data/stress-test-0/test-4-polygons-abcd.off", { 1, 2 }, results);
-  results[0] = { 54, 76, 16 };
-  results[1] = { 64, 102, 24 };
+  run_test<Kernel, IntersectionKernel>("data/stress-test-0/test-4-polygons-abcd.off", { 1, 2 }, results);
+  results[0] = { 50, 68, 14 };
+  results[1] = { 56, 82, 18 };
   results[2] = { 67, 109, 26 };
-  run_test<Kernel>("data/stress-test-0/test-6-polygons.off", { 1, 2, 3 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-0/test-6-polygons.off", { 1, 2, 3 }, results);
 
-  // Stress tests 1.
+ // Stress tests 1.
 
   results[0] = { 14, 13, 2 };
-  run_test<Kernel>("data/stress-test-1/test-1-rnd-polygons-1-4.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-1/test-1-rnd-polygons-1-4.off", { 1 }, results);
   results[0] = { 14, 13, 2 };
-  run_test<Kernel>("data/stress-test-1/test-2-rnd-polygons-1-4.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-1/test-2-rnd-polygons-1-4.off", { 1 }, results);
   results[0] = { 14, 13, 2 };
-  run_test<Kernel>("data/stress-test-1/test-3-rnd-polygons-1-4.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-1/test-3-rnd-polygons-1-4.off", { 1 }, results);
   results[0] = { 14, 13, 2 };
-  run_test<Kernel>("data/stress-test-1/test-4-rnd-polygons-1-4.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-1/test-4-rnd-polygons-1-4.off", { 1 }, results);
   results[0] = { 18, 18, 3 };
   results[1] = { 20, 22, 4 };
-  run_test<Kernel>("data/stress-test-1/test-5-rnd-polygons-2-4.off", { 1, 2 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-1/test-5-rnd-polygons-2-4.off", { 1, 2 }, results);
   results[0] = { 18, 18, 3 };
   results[1] = { 19, 21, 4 };
-  run_test<Kernel>("data/stress-test-1/test-6-rnd-polygons-2-4.off", { 1, 2 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-1/test-6-rnd-polygons-2-4.off", { 1, 2 }, results);
   results[0] = { 18, 18, 3 };
-  run_test<Kernel>("data/stress-test-1/test-7-rnd-polygons-2-4.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-1/test-7-rnd-polygons-2-4.off", { 1 }, results);
 
   // Stress tests 2.
 
   results[0] = { 14, 13, 2 };
-  run_test<Kernel>("data/stress-test-2/test-1-rnd-polygons-1-4.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-2/test-1-rnd-polygons-1-4.off", { 1 }, results);
   results[0] = { 14, 13, 2 };
-  run_test<Kernel>("data/stress-test-2/test-2-rnd-polygons-1-4.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-2/test-2-rnd-polygons-1-4.off", { 1 }, results);
   results[0] = { 14, 13, 2 };
-  run_test<Kernel>("data/stress-test-2/test-3-rnd-polygons-1-4.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-2/test-3-rnd-polygons-1-4.off", { 1 }, results);
   results[0] = { 14, 13, 2 };
-  run_test<Kernel>("data/stress-test-2/test-4-rnd-polygons-1-3.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-2/test-4-rnd-polygons-1-3.off", { 1 }, results);
   results[0] = { 17, 17, 3 };
   results[1] = { 19, 21, 4 };
-  run_test<Kernel>("data/stress-test-2/test-5-rnd-polygons-2-4.off", { 1, 2 }, results);
-  results[0] = { 22, 23, 4 };
-  run_test<Kernel>("data/stress-test-2/test-6-rnd-polygons-3-4.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-2/test-5-rnd-polygons-2-4.off", { 1, 2 }, results);
+  results[0] = { 24, 27, 5 };
+  run_test<Kernel, IntersectionKernel>("data/stress-test-2/test-6-rnd-polygons-3-4.off", { 1 }, results);
 
   // Stress tests 3.
 
   results[0] = { 18, 18, 3 };
   results[1] = { 20, 22, 4 };
-  run_test<Kernel>("data/stress-test-3/test-1-rnd-polygons-2-3.off", { 1, 2 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-3/test-1-rnd-polygons-2-3.off", { 1, 2 }, results);
   results[0] = { 17, 17, 3 };
-  run_test<Kernel>("data/stress-test-3/test-2-rnd-polygons-2-3.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-3/test-2-rnd-polygons-2-3.off", { 1 }, results);
   results[0] = { 18, 18, 3 };
   results[1] = { 19, 21, 4 };
-  run_test<Kernel>("data/stress-test-3/test-3-rnd-polygons-2-3.off", { 1, 2 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-3/test-3-rnd-polygons-2-3.off", { 1, 2 }, results);
   results[0] = { 17, 17, 3 };
   results[1] = { 19, 21, 4 };
-  run_test<Kernel>("data/stress-test-3/test-4-rnd-polygons-2-4.off", { 1, 2 }, results);
-  //results[0] = { 12, 11, 2 };
-  //run_test<Kernel>("data/stress-test-3/test-5-rnd-polygons-1-3.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-3/test-4-rnd-polygons-2-4.off", { 1, 2 }, results);
+  results[0] = { 14, 13, 2 };
+  run_test<Kernel, IntersectionKernel>("data/stress-test-3/test-5-rnd-polygons-1-3.off", { 1 }, results);
   results[0] = { 18, 18, 3 };
   results[1] = { 19, 21, 4 };
-  run_test<Kernel>("data/stress-test-3/test-6-rnd-polygons-2-3.off", { 1, 2 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-3/test-6-rnd-polygons-2-3.off", { 1, 2 }, results);
   results[0] = { 21, 21, 3 };
   results[1] = { 22, 24, 4 };
-  run_test<Kernel>("data/stress-test-3/test-7-rnd-polygons-2-4.off", { 1, 2 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-3/test-7-rnd-polygons-2-4.off", { 1, 2 }, results);
   results[0] = { 17, 17, 3 };
   results[1] = { 18, 20, 4 };
-  run_test<Kernel>("data/stress-test-3/test-8-rnd-polygons-2-10.off", { 1, 2 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-3/test-8-rnd-polygons-2-10.off", { 1, 2 }, results);
   results[0] = { 31, 37, 7 };
   results[1] = { 34, 46, 10 };
   results[2] = { 39, 57, 13 };
-  run_test<Kernel>("data/stress-test-3/test-9-rnd-polygons-4-4.off", { 1, 2, 3 }, results);
-  results[0] = { 55, 84, 19 };
-  run_test<Kernel>("data/stress-test-3/test-10-rnd-polygons-5-4.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-3/test-9-rnd-polygons-4-4.off", { 1, 2, 3 }, results);
+  results[0] = { 51, 72, 15 };
+  run_test<Kernel, IntersectionKernel>("data/stress-test-3/test-10-rnd-polygons-5-4.off", { 1 }, results);
 
   // Stress tests 4.
 
   results[0] = { 17, 17, 3 };
-  run_test<Kernel>("data/stress-test-4/test-1-rnd-polygons-2-6.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-4/test-1-rnd-polygons-2-6.off", { 1 }, results);
   results[0] = { 27, 32, 6 };
   results[1] = { 29, 38, 8 };
-  run_test<Kernel>("data/stress-test-4/test-2-rnd-polygons-3-8.off", { 1, 2 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-4/test-2-rnd-polygons-3-8.off", { 1, 2 }, results);
   results[0] = { 32, 38, 7 };
   results[1] = { 35, 45, 9 };
   results[2] = { 37, 51, 11 };
-  run_test<Kernel>("data/stress-test-4/test-3-rnd-polygons-4-4.off", { 1, 2, 3 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-4/test-3-rnd-polygons-4-4.off", { 1, 2, 3 }, results);
   results[0] = { 25, 27, 5 };
-  results[1] = { 33, 41, 8 };
+  results[1] = { 30, 36, 7 };
   results[2] = { 37, 53, 12 };
-  run_test<Kernel>("data/stress-test-4/test-4-rnd-polygons-4-6.off", { 1, 2, 3 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-4/test-4-rnd-polygons-4-6.off", { 1, 2, 3 }, results);
   results[0] = { 61, 91, 20 };
   results[1] = { 73, 121, 29 };
   results[2] = { 83, 145, 36 };
-  run_test<Kernel>("data/stress-test-4/test-5-rnd-polygons-6-4.off", { 1, 2, 3 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-4/test-5-rnd-polygons-6-4.off", { 1, 2, 3 }, results);
 
-  results[0] = { 45, 62, 13 };
+  results[0] = { 43, 58, 12 };
   results[1] = { 50, 75, 17 };
-  run_test<Kernel>("data/stress-test-4/test-6-rnd-polygons-5-6.off", { 1, 2 }, results);
-  results[0] = { 64, 97, 22 };
+  run_test<Kernel, IntersectionKernel>("data/stress-test-4/test-6-rnd-polygons-5-6.off", { 1, 2 }, results);
+  results[0] = { 63, 94, 21 };
   results[1] = { 84, 141, 34 };
-  results[2] = { 88, 151, 37 };
-  run_test<Kernel>("data/stress-test-4/test-7-rnd-polygons-7-6.off", { 1, 2, 3 }, results);
+  results[2] = { 90, 157, 39 };
+  run_test<Kernel, IntersectionKernel>("data/stress-test-4/test-7-rnd-polygons-7-6.off", { 1, 2, 3 }, results);
   results[0] = { 56, 77, 16 };
   results[1] = { 68, 107, 25 };
   results[2] = { 69, 110, 26 };
-  run_test<Kernel>("data/stress-test-4/test-8-rnd-polygons-7-8.off", { 1, 2, 3 }, results);
-  results[0] = { 172, 304, 74 };
-  results[1] = { 192, 366, 95 };
-  results[2] = { 198, 382, 100};
-  run_test<Kernel>("data/stress-test-4/test-9-rnd-polygons-12-4.off", { 1, 2, 3 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-4/test-8-rnd-polygons-7-8.off", { 1, 2, 3 }, results);
+  results[0] = { 173, 305, 74 };
+  results[1] = { 194, 370, 96 };
+  results[2] = { 207, 407, 108 };
+  run_test<Kernel, IntersectionKernel>("data/stress-test-4/test-9-rnd-polygons-12-4.off", { 1, 2, 3 }, results);
 
   // Stress tests 5.
 
-  results[0] = { 202, 351, 84 };
-  results[1] = { 232, 427, 107 };
-  results[2] = { 284, 558, 146 };
-  run_test<Kernel>("data/stress-test-5/test-1-rnd-polygons-15-6.off", { 1, 2, 3 }, results);
-  results[0] = { 58, 89, 20 };
-  results[1] = { 63, 102, 24 };
-  results[2] = { 63, 106, 26 };
-  run_test<Kernel>("data/stress-test-5/test-2-rnd-polygons-20-4.off", { 1, 2, 3 }, results);
+  results[0] = { 185, 308, 71 };
+  results[1] = { 223, 406, 101 };
+  results[2] = { 277, 548, 145 };
+  run_test<Kernel, IntersectionKernel>("data/stress-test-5/test-1-rnd-polygons-15-6.off", { 1, 2, 3 }, results);
+
+  results[0] = { 50, 71, 15 };
+  results[1] = { 56, 85, 19 };
+  results[2] = { 63, 102, 24 };
+  run_test<Kernel, IntersectionKernel>("data/stress-test-5/test-2-rnd-polygons-20-4.off", { 1, 2, 3 }, results);
 
   // Real data tests.
 
-  results[0] = { 91, 143, 33 };
-  results[1] = { 109, 187, 46 };
+  results[0] = { 94, 150, 35 };
+  results[1] = { 111, 191, 47 };
   results[2] = { 127, 233, 60 };
-  run_test<Kernel>("data/real-data-test/test-10-polygons.off", { 1, 2, 3 }, results);
+  run_test<Kernel, IntersectionKernel>("data/real-data-test/test-10-polygons.off", { 1, 2, 3 }, results);
 
-  results[0] = { 1006, 2067, 552 };
-  results[1] = { 973, 1984, 527 };
-  results[2] = { 1186, 2560, 708 };
-  run_test<Kernel>("data/real-data-test/test-40-polygons.ply", { 1, 2, 3 }, results);
+  results[0] = { 206, 385, 99 };
+  results[1] = { 237, 462, 122 };
+  results[2] = { 260, 529, 144 };
+  run_test<Kernel, IntersectionKernel>("data/real-data-test/test-15-polygons.off", { 1, 2, 3 }, results);
+
+  results[0] = { 1156, 2466, 677 };
+  results[1] = { 1131, 2387, 650 };
+  results[2] = { 1395, 3115, 882 };
+  run_test<Kernel, IntersectionKernel>("data/real-data-test/test-40-polygons.ply", { 1, 2, 3 }, results);
 
   const auto kernel_name = boost::typeindex::type_id<Kernel>().pretty_name();
+  if (different != 0) {
+    std::cout << std::endl << kernel_name << different << " TESTS differ from typical values!" << std::endl << std::endl;
+  }
   std::cout << std::endl << kernel_name << " TESTS SUCCESS!" << std::endl << std::endl;
 }
 
-#include <CGAL/intersections.h>
-
 int main(const int /* argc */, const char** /* argv */) {
-  // run_all_tests<SCF>();
-  // run_all_tests<SCD>();
-  //run_all_tests<EPECK>();
-
-  // Passes all tests except for those when
-  // intersections lead to accumulated errors.
-  //build();
-  run_all_tests<EPICK>();
+  run_all_tests<SCD, EPECK>();
+  run_all_tests<EPICK, EPECK>();
+  //run_all_tests<EPECK, EPECK>();
+  //run_all_tests<GMPQ, GMPQ>();
   return EXIT_SUCCESS;
 }