diff --git a/Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/CMakeLists.txt b/Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/CMakeLists.txt
index b4737eceb39..246a2b3ca3b 100644
--- a/Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/CMakeLists.txt
+++ b/Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/CMakeLists.txt
@@ -19,11 +19,11 @@ if(CGAL_FOUND)
   if(Boost_FOUND)
     message(STATUS "Found Boost")
 
-    # find_package(Eigen3 3.1.0 REQUIRED)
-    # if(Eigen3_FOUND)
-    #   message(STATUS "Found Eigen")
+    find_package(Eigen3 3.1.0 REQUIRED)
+    if(Eigen3_FOUND)
+      message(STATUS "Found Eigen")
 
-    #   include(CGAL_Eigen_support)
+      include(CGAL_Eigen_support)
 
       set(targets
       # kinetic_2d_example
@@ -37,13 +37,13 @@ if(CGAL_FOUND)
       foreach(target ${targets})
         create_single_source_cgal_program("${target}.cpp")
         if(TARGET ${target})
-          target_link_libraries(${target} PUBLIC ${project_linked_libraries})
+          target_link_libraries(${target} PUBLIC ${project_linked_libraries} CGAL::Eigen_support)
           target_compile_definitions(${target} PUBLIC ${project_compilation_definitions})
         endif()
       endforeach()
-    # else()
-    #   message(ERROR "This program requires the Eigen library, and will not be compiled.")
-    # endif()
+    else()
+      message(ERROR "This program requires the Eigen library, and will not be compiled.")
+    endif()
   else()
     message(ERROR "This program requires the Boost library, and will not be compiled.")
   endif()
diff --git a/Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/data/edge-case-test/test-flat-bbox.off b/Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/data/edge-case-test/test-flat-bbox.off
new file mode 100644
index 00000000000..f36cb73dce9
--- /dev/null
+++ b/Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/data/edge-case-test/test-flat-bbox.off
@@ -0,0 +1,12 @@
+OFF
+8 2 0
+0.0 0.0 0.0
+1.0 0.0 0.0
+1.0 1.0 0.0
+0.0 1.0 0.0
+2.0 0.0 0.0
+3.0 0.0 0.0
+3.0 1.0 0.0
+2.0 1.0 0.0
+4 0 1 2 3
+4 4 5 6 7
diff --git a/Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/kinetic_precomputed_shapes_example.cpp b/Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/kinetic_precomputed_shapes_example.cpp
index 00ab562fb27..1a083c6419a 100644
--- a/Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/kinetic_precomputed_shapes_example.cpp
+++ b/Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/kinetic_precomputed_shapes_example.cpp
@@ -70,14 +70,15 @@ int main(const int argc, const char** argv) {
   const unsigned int k = (argc > 2 ? std::atoi(argv[2]) : 1);
   std::cout << "* number of intersections k: " << k << std::endl;
 
-  const unsigned int n = 0; // number of subdivisions per bbox side
+  const unsigned int n = 0;
+  std::cout << "* number of subdivisions per bbox side: " << n << std::endl;
   const unsigned int num_blocks = std::pow(n + 1, 3);
   std::cout << "* number of blocks: " << num_blocks << std::endl;
 
   const double enlarge_bbox_ratio = 1.1;
   std::cout << "* enlarge bbox ratio: " << enlarge_bbox_ratio << std::endl;
 
-  const bool reorient = true;
+  const bool reorient = false;
   std::cout << "* reorient: " << (reorient ? "true" : "false") << std::endl;
 
   // Algorithm.
diff --git a/Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/kinetic_random_shapes_example.cpp b/Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/kinetic_random_shapes_example.cpp
index 57edc52e5b1..0f93f69f3f4 100644
--- a/Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/kinetic_random_shapes_example.cpp
+++ b/Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/kinetic_random_shapes_example.cpp
@@ -380,11 +380,11 @@ int main(const int argc, const char** argv) {
 
   std::cout << std::endl;
   std::cout << "--- INPUT STATS: " << std::endl;
-  std::cout << "* input kernel: " << boost::typeindex::type_id<EPICK>().pretty_name() << std::endl;
-  std::cout << "* polygon kernel: " << boost::typeindex::type_id<Kernel>().pretty_name() << std::endl;
-  std::cout << "* expected number of polygons: " << n << std::endl;
-  std::cout << "* generated number of polygons: " << rnd_polygons.size() << std::endl;
-  std::cout << "* number of vertices in a polygon: " << p << std::endl;
+  std::cout << "* input kernel: "                    << boost::typeindex::type_id<EPICK>().pretty_name()  << std::endl;
+  std::cout << "* polygon kernel: "                  << boost::typeindex::type_id<Kernel>().pretty_name() << std::endl;
+  std::cout << "* expected number of polygons: "     << n                                                 << std::endl;
+  std::cout << "* generated number of polygons: "    << rnd_polygons.size()                               << std::endl;
+  std::cout << "* number of vertices in a polygon: " << p                                                 << std::endl;
   // exit(EXIT_SUCCESS);
 
   IPolygon_3 input_polygon;
@@ -403,7 +403,7 @@ int main(const int argc, const char** argv) {
   assert(input_polygons.size() == rnd_polygons.size());
 
   // Algorithm.
-  KSR ksr(true, true);
+  KSR ksr(false, false);
   const IPolygon_3_map polygon_map;
   const unsigned int k = (argc > 3 ? std::atoi(argv[3]) : 1);
   std::cout << "* input k: " << k << std::endl;
diff --git a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Data_structure.h b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Data_structure.h
index 640fe394e41..afa25514784 100644
--- a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Data_structure.h
+++ b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Data_structure.h
@@ -2422,11 +2422,11 @@ public:
         pvertex, prev, next, crossed[i], future_points[i], future_directions[i]);
       if (is_parallel) {
         if (is_intersecting_iedge(min_time, max_time, prev, crossed[i])) {
-          CGAL_assertion_msg(i == crossed.size() - 1, "TODO: FRONT, CAN WE HAVE OTHER I HERE?");
+          CGAL_assertion_msg(i == crossed.size() - 1, "TODO: OPEN, PREV, CAN WE HAVE OTHER I HERE?");
           prev_iedge = crossed[i];
         }
         if (is_intersecting_iedge(min_time, max_time, next, crossed[i])) {
-          CGAL_assertion_msg(i == 0, "TODO: FRONT, CAN WE HAVE OTHER I HERE?");
+          CGAL_assertion_msg(i == 0, "TODO: OPEN, NEXT, CAN WE HAVE OTHER I HERE?");
           next_iedge = crossed[i];
         }
       }
@@ -2654,7 +2654,6 @@ public:
         CGAL_assertion(this->k(pface_front) == 1);
         CGAL_assertion(this->k(pface_back)  == 1);
       }
-      // CGAL_assertion_msg(false, "TODO: FRONT AND BACK, FINISH THIS CASE!");
 
     } else if ((!is_occupied_edge_front && !is_occupied_edge_back)) {
 
@@ -2662,9 +2661,7 @@ public:
       CGAL_assertion(this->k(pface_front) == this->k(pface_back));
       add_new_pfaces(this->k(pface_front), pvertex, crossed,
       future_points, future_directions, new_pvertices);
-
       if (m_verbose) std::cout << "- continue !front && !back" << std::endl;
-      // CGAL_assertion_msg(false, "TODO: !FRONT AND !BACK, FINISH THIS CASE!");
 
     } else if (is_occupied_edge_front || is_occupied_edge_back) {
 
@@ -2681,9 +2678,7 @@ public:
       } else {
         CGAL_assertion_msg(false, "ERROR: WRONG OPEN CASE!");
       }
-
       if (m_verbose) std::cout << "- continue front || back" << std::endl;
-      // CGAL_assertion_msg(false, "TODO: FRONT OR BACK, FINISH THIS CASE!");
 
     } else {
       CGAL_assertion_msg(false, "TODO: ADD NEW OPEN CASE! DO NOT FORGET TO UPDATE K!");
@@ -3608,7 +3603,7 @@ public:
       auto vec2 = Vector_3(centroid, other_centroid);
       vec2 = KSR::normalize(vec2);
 
-      const FT d = FT(1) / FT(100000); // TODO: CAN WE AVOID THIS VALUE?
+      const FT d = KSR::tolerance<FT>(); // TODO: CAN WE AVOID THIS VALUE?
       const FT dot_product = vec1 * vec2;
 
       if (dot_product < FT(0)) {
diff --git a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Event.h b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Event.h
index f67e5da20e7..8502ae5383b 100644
--- a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Event.h
+++ b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Event.h
@@ -28,7 +28,7 @@
 namespace CGAL {
 namespace KSR_3 {
 
-// TODO: Can we avoid forward declaration?
+// TODO: CAN WE AVOID FORWARD DECLARATION?
 template<typename Data_structure>
 class Event_queue;
 
diff --git a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Event_queue.h b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Event_queue.h
index c1f2fc0a9f0..655484c8d3a 100644
--- a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Event_queue.h
+++ b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Event_queue.h
@@ -146,7 +146,7 @@ public:
     }
     queue_by_pvertex_idx().erase(pv.first, pv.second);
 
-    // Erase by pother. TODO: Why is pother here?
+    // Erase by pother.
     const auto po = queue_by_pother_idx().equal_range(pvertex);
     const auto po_range = CGAL::make_range(po);
 
diff --git a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Experimental.h b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Experimental.h
index 4a90fd1f949..433936b25d9 100644
--- a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Experimental.h
+++ b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Experimental.h
@@ -560,7 +560,7 @@ class Experimental {
         // continue ..
         // std::cout << "crossed size: " << crossed.size() << std::endl;
         // std::cout << "all crossed size: " << all_crossed.size() << std::endl;
-        // CGAL_assertion_msg(false, "TODO: BACK CROSSED > LIMIT!");
+        CGAL_assertion_msg(false, "TODO: BACK CROSSED > LIMIT!");
       }
     }
     else if (front_constrained) // Border case
@@ -893,7 +893,7 @@ class Experimental {
         if (is_ok) {
         if (num_extra_faces < 3) {
 
-          // CGAL_assertion_msg(false, "TODO: FRONT, CROSSED < LIMIT, 1 or 2 FACES!");
+          CGAL_assertion_msg(false, "TODO: FRONT, CROSSED < LIMIT, 1 or 2 FACES!");
           CGAL_assertion(future_points.size() == asize);
           CGAL_assertion(future_directions.size() == asize);
 
@@ -943,7 +943,7 @@ class Experimental {
               connect(propagated, all_crossed[i]);
               crossed.push_back(all_crossed[i]); // remove events from this one
 
-              // CGAL_assertion_msg(false, "TODO: FRONT, NULL PROPAGATED CASE!");
+              CGAL_assertion_msg(false, "TODO: FRONT, NULL PROPAGATED CASE!");
 
             } else {
 
@@ -1116,7 +1116,7 @@ class Experimental {
         // continue ..
         // std::cout << "crossed size: " << crossed.size() << std::endl;
         // std::cout << "all crossed size: " << all_crossed.size() << std::endl;
-        // CGAL_assertion_msg(false, "TODO: FRONT, CROSSED > LIMIT!");
+        CGAL_assertion_msg(false, "TODO: FRONT, CROSSED > LIMIT!");
       }
     }
     else // Open case
diff --git a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Initializer.h b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Initializer.h
index d34de855e06..fb00b18e9c2 100644
--- a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Initializer.h
+++ b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Initializer.h
@@ -24,7 +24,8 @@
 // #include <CGAL/license/Kinetic_shape_reconstruction.h>
 
 // CGAL includes.
-// #include <CGAL/optimal_bounding_box.h>
+#include <CGAL/optimal_bounding_box.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 
 // Internal includes.
 #include <CGAL/KSR/utils.h>
@@ -51,10 +52,13 @@ private:
   using Data_structure   = KSR_3::Data_structure<Kernel>;
   using Polygon_splitter = KSR_3::Polygon_splitter<Data_structure, Kernel>;
 
-  using IVertex = typename Data_structure::IVertex;
+  using IVertex  = typename Data_structure::IVertex;
+  using IK       = CGAL::Exact_predicates_inexact_constructions_kernel;
+  using IFT      = typename IK::FT;
+  using IPoint_3 = typename IK::Point_3;
 
-  using Bbox_3 = CGAL::Bbox_3;
-  // using OBB_traits = CGAL::Oriented_bounding_box_traits_3<Kernel>;
+  using Bbox_3     = CGAL::Bbox_3;
+  using OBB_traits = CGAL::Oriented_bounding_box_traits_3<IK>;
 
 public:
   Initializer(
@@ -155,10 +159,11 @@ private:
     std::array<Point_3, 8>& bbox,
     FT& time_step) const {
 
-    if (reorient)
+    if (reorient) {
       initialize_optimal_box(input_range, polygon_map, bbox);
-    else
+    } else {
       initialize_axis_aligned_box(input_range, polygon_map, bbox);
+    }
 
     CGAL_assertion(bbox.size() == 8);
     time_step  = KSR::distance(bbox.front(), bbox.back());
@@ -194,24 +199,55 @@ private:
     }
 
     // Set points.
-    std::vector<Point_3> points;
-    points.reserve(num_points);
+    std::vector<IPoint_3> ipoints;
+    ipoints.reserve(num_points);
     for (const auto& item : input_range) {
       const auto& polygon = get(polygon_map, item);
-      for (const auto& p : polygon) {
-        const Point_3 point(p.x(), p.y(), p.z());
-        points.push_back(point);
+      for (const auto& point : polygon) {
+        const IPoint_3 ipoint(
+          static_cast<IFT>(CGAL::to_double(point.x())),
+          static_cast<IFT>(CGAL::to_double(point.y())),
+          static_cast<IFT>(CGAL::to_double(point.z())));
+        ipoints.push_back(ipoint);
       }
     }
 
     // Compute optimal bbox.
     // The order of faces corresponds to the standard order from here:
     // https://doc.cgal.org/latest/BGL/group__PkgBGLHelperFct.html#gad9df350e98780f0c213046d8a257358e
-    // const OBB_traits obb_traits;
-    // CGAL::oriented_bounding_box(
-    //   points, bbox,
-    //   CGAL::parameters::use_convex_hull(true).
-    //   geom_traits(obb_traits));
+    const OBB_traits obb_traits;
+    std::array<IPoint_3, 8> ibbox;
+    CGAL::oriented_bounding_box(
+      ipoints, ibbox,
+      CGAL::parameters::use_convex_hull(true).
+      geom_traits(obb_traits));
+
+    for (std::size_t i = 0; i < 8; ++i) {
+      const auto& ipoint = ibbox[i];
+      const Point_3 point(
+        static_cast<FT>(ipoint.x()),
+        static_cast<FT>(ipoint.y()),
+        static_cast<FT>(ipoint.z()));
+      bbox[i] = point;
+    }
+
+    const FT sq_length = CGAL::squared_distance(bbox[0], bbox[1]);
+    const FT sq_width  = CGAL::squared_distance(bbox[0], bbox[3]);
+    const FT sq_height = CGAL::squared_distance(bbox[0], bbox[5]);
+    CGAL_assertion(sq_length >= FT(0));
+    CGAL_assertion(sq_width  >= FT(0));
+    CGAL_assertion(sq_height >= FT(0));
+    const FT tol = KSR::tolerance<FT>();
+    if (sq_length < tol || sq_width < tol || sq_height < tol) {
+      if (m_verbose) {
+        std::cout << "* warning: optimal bounding box is flat, reverting ..." << std::endl;
+      }
+      initialize_axis_aligned_box(input_range, polygon_map, bbox);
+    } else {
+      if (m_verbose) {
+        std::cout << "* using optimal bounding box" << std::endl;
+      }
+    }
   }
 
   template<
@@ -239,16 +275,34 @@ private:
       Point_3(box.xmin(), box.ymin(), box.zmax()),
       Point_3(box.xmax(), box.ymin(), box.zmax()),
       Point_3(box.xmax(), box.ymax(), box.zmax()) };
+
+    const FT sq_length = CGAL::squared_distance(bbox[0], bbox[1]);
+    const FT sq_width  = CGAL::squared_distance(bbox[0], bbox[3]);
+    const FT sq_height = CGAL::squared_distance(bbox[0], bbox[5]);
+    CGAL_assertion(sq_length >= FT(0));
+    CGAL_assertion(sq_width  >= FT(0));
+    CGAL_assertion(sq_height >= FT(0));
+    const FT tol = KSR::tolerance<FT>();
+    if (sq_length < tol || sq_width < tol || sq_height < tol) {
+      CGAL_assertion_msg(false, "TODO: HANDLE FLAT AXIS ALIGNED BBOX!");
+    } else {
+      if (m_verbose) {
+        std::cout << "* using axis aligned bounding box" << std::endl;
+      }
+    }
   }
 
   void enlarge_bounding_box(
     const FT enlarge_bbox_ratio,
     std::array<Point_3, 8>& bbox) const {
 
-    CGAL_assertion_msg(
-      enlarge_bbox_ratio > FT(1), "TODO: HANDLE THE CASE ENLARGE_BBOX_RATIO = FT(1)");
+    FT enlarge_ratio = enlarge_bbox_ratio;
+    if (enlarge_bbox_ratio == FT(1)) {
+      enlarge_ratio += KSR::tolerance<FT>();
+    }
+
     const auto a = CGAL::centroid(bbox.begin(), bbox.end());
-    Transform_3 scale(CGAL::Scaling(), enlarge_bbox_ratio);
+    Transform_3 scale(CGAL::Scaling(), enlarge_ratio);
     for (auto& point : bbox)
       point = scale.transform(point);
 
diff --git a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Polygon_splitter.h b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Polygon_splitter.h
index b4cc26ff52f..b2b73d10c6d 100644
--- a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Polygon_splitter.h
+++ b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Polygon_splitter.h
@@ -25,7 +25,6 @@
 
 // CGAL includes.
 #include <CGAL/convex_hull_2.h>
-#include <CGAL/Regular_triangulation_2.h>
 #include <CGAL/Triangulation_face_base_with_info_2.h>
 #include <CGAL/Constrained_Delaunay_triangulation_2.h>
 #include <CGAL/Constrained_triangulation_plus_2.h>
@@ -98,11 +97,6 @@ private:
   using Face_index   = typename Mesh_3::Face_index;
   using Uchar_map    = typename Mesh_3::template Property_map<Face_index, unsigned char>;
 
-  using Regular_triangulation = CGAL::Regular_triangulation_2<Kernel>;
-  using Weighted_point        = typename Regular_triangulation::Weighted_point;
-  using RVertex_handle        = typename Regular_triangulation::Vertex_handle;
-  using REdge                 = typename Regular_triangulation::Edge;
-
   enum struct Merge_type { CONVEX_HULL = 0, RECTANGLE = 1 };
 
   Data_structure& m_data;
@@ -147,9 +141,6 @@ public:
     if (original_faces.size() > 1) {
       CGAL_assertion_msg(false,
       "ERROR: WE CANNOT HAVE MULTIPLE COPLANAR PFACES!");
-      // merge_coplanar_pfaces(
-      //   support_plane_idx, original_input, original_faces);
-      // add_bisectors(original_faces, support_plane_idx);
     }
   }
 
@@ -180,9 +171,7 @@ private:
         m_data.to_3d(support_plane_idx, q) << std::endl;
       }
     }
-
     add_merged_pface(support_plane_idx, merged);
-    // CGAL_assertion_msg(false, "TODO: MERGE COPLANAR PFACES!");
   }
 
   void collect_pface_points(
@@ -234,14 +223,14 @@ private:
     const KSR::size_t support_plane_idx,
     const std::vector<Point_2>& merged) const {
 
-    std::vector<Weighted_point> wps;
-    add_weighted_bbox(support_plane_idx, wps);
-    CGAL_assertion(wps.size() == 4);
+    std::vector<Point_2> bbox;
+    create_bbox(support_plane_idx, bbox);
+    CGAL_assertion(bbox.size() == 4);
 
     for (std::size_t i = 0; i < 4; ++i) {
       const std::size_t ip = (i + 1) % 4;
-      const auto& pi = wps[i].point();
-      const auto& qi = wps[ip].point();
+      const auto& pi = bbox[i];
+      const auto& qi = bbox[ip];
       const Segment_2 edge(pi, qi);
 
       for (std::size_t j = 0; j < merged.size(); ++j) {
@@ -295,30 +284,25 @@ private:
     original_faces.clear();
     original_input.clear();
 
-    // std::vector<Vertex_handle> vhs; // TODO: Why we cannot use vhs here?
-    std::vector<TRI> triangulations;
+    // std::vector<Vertex_handle> vhs;
     std::vector<Point_2> original_face;
-
     const auto all_pfaces = m_data.pfaces(support_plane_idx);
     for (const auto pface : all_pfaces) {
       const auto pvertices = m_data.pvertices_of_pface(pface);
 
       // vhs.clear();
-      // TRI triangulation; // TODO: remove triangulations here!
       original_face.clear();
       for (const auto pvertex : pvertices) {
         CGAL_assertion(vhs_map.find(pvertex) != vhs_map.end());
         const auto vh = vhs_map.at(pvertex);
         original_face.push_back(vh->point());
-        // triangulation.insert(vh->point());
         // vhs.push_back(vh);
       }
 
-      // triangulations.push_back(triangulation);
       original_faces.push_back(original_face);
       original_input.push_back(m_data.input(pface));
 
-      // TODO: Why we cannot use vhs directly here? That should be more precise!
+      // TODO: WHY WE CANNOT USE VHS DIRECTLY HERE? THAT SHOULD BE MORE PRECISE!
       // vhs.push_back(vhs.front());
       // const auto cid = m_cdt.insert_constraint(vhs.begin(), vhs.end());
       original_face.push_back(original_face.front());
@@ -343,21 +327,8 @@ private:
 
     // Finally, add original labels to the cdt.
     if (all_pfaces.size() > 1) {
-      CGAL_assertion_msg(false, "ERROR: WE CANNOT HAVE MULTIPLE COPLANAR PFACES!");
-      // for (auto fit = m_cdt.finite_faces_begin(); fit != m_cdt.finite_faces_end(); ++fit) {
-      //   const auto centroid = CGAL::centroid(
-      //     fit->vertex(0)->point(),
-      //     fit->vertex(1)->point(),
-      //     fit->vertex(2)->point());
-      //   for (KSR::size_t i = 0; i < triangulations.size(); ++i) {
-      //     const auto fh = triangulations[i].locate(centroid);
-      //     if (fh == nullptr || triangulations[i].is_infinite(fh)) {
-      //       continue;
-      //     }
-      //     fit->info().input = i;
-      //     break;
-      //   }
-      // }
+      CGAL_assertion_msg(false,
+      "ERROR: WE CANNOT HAVE MULTIPLE COPLANAR PFACES!");
     }
   }
 
@@ -525,16 +496,10 @@ private:
       if (original_faces.size() != 1) {
         CGAL_assertion_msg(original_faces.size() <= 1,
         "ERROR: WE CANNOT HAVE MULTIPLE COPLANAR PFACES!");
-        // dump_original_faces(support_plane_idx, original_faces, "original-faces");
-        // dump_current_pface(pface, "current-pface-" + m_data.str(pface));
-        // locate_pface_among_coplanar_pfaces(pface);
       } else {
         m_data.input(pface) = original_input[0];
       }
     }
-    // if (original_faces.size() > 1) {
-    //   CGAL_assertion_msg(false, "TODO: DEBUG THIS SUPPORT PLANE!");
-    // }
   }
 
   void reconnect_pvertices_to_ivertices() {
@@ -695,159 +660,9 @@ private:
     return is_okay;
   }
 
-  void locate_pface_among_coplanar_pfaces(
-    const PFace& pface) {
-
-    // TODO: Change fh->info().input to vector!
-    // std::cout << "locating " << m_data.str(pface) << std::endl;
-    const auto centroid = m_data.centroid_of_pface(pface);
-    const auto fh = m_cdt.locate(m_data.to_2d(pface.first, centroid));
-    CGAL_assertion(fh != nullptr && !m_cdt.is_infinite(fh));
-    CGAL_assertion(fh->info().input != KSR::uninitialized());
-    m_data.input(pface) = fh->info().input;
-    // std::cout << "found original input: " << m_data.input(pface) << std::endl;
-  }
-
-  void merge_coplanar_pfaces(
+  void create_bbox(
     const KSR::size_t support_plane_idx,
-    const std::vector< std::vector<KSR::size_t> >& original_input,
-    const std::vector< std::vector<Point_2> >& original_faces) {
-
-    const bool is_debug = false;
-    CGAL_assertion(support_plane_idx >= 6);
-
-    if (is_debug) {
-      std::cout << std::endl << "support plane idx: " << support_plane_idx << std::endl;
-      std::cout << "dumping support plane ... ";
-      dump(true, 0, support_plane_idx, "support-plane-" +
-      std::to_string(support_plane_idx) + ".ply");
-      std::cout << "done" << std::endl;
-    }
-
-    if (original_faces.size() > 2) {
-      CGAL_assertion_msg(false, "TODO: MERGE > 2 COPLANAR PFACES!");
-      // This code should be very similar to the one with 2 coplanar pfaces!
-    }
-
-    const auto all_pfaces = m_data.pfaces(support_plane_idx);
-    std::vector<PFace> pfaces;
-    pfaces.reserve(all_pfaces.size());
-    for (const auto pface : all_pfaces) {
-      pfaces.push_back(pface);
-    }
-    CGAL_assertion(pfaces.size() >= 2);
-    CGAL_assertion(pfaces.size() == all_pfaces.size());
-    if (is_debug) std::cout << "num pfaces: " << pfaces.size() << std::endl;
-
-    std::vector< std::pair<PFace, PFace> > to_be_merged;
-    const auto& iedges = m_data.iedges(support_plane_idx);
-    for (std::size_t i = 0; i < pfaces.size(); ++i) {
-      const auto& pface_i = pfaces[i];
-      const auto centroid_i = m_data.centroid_of_pface(pface_i);
-      for (std::size_t j = i + 1; j < pfaces.size(); ++j) {
-        const auto& pface_j = pfaces[j];
-        const auto centroid_j = m_data.centroid_of_pface(pface_j);
-
-        const Segment_2 query_segment(
-          m_data.to_2d(support_plane_idx, centroid_i),
-          m_data.to_2d(support_plane_idx, centroid_j));
-
-        bool found_intersection = false;
-        for (const auto& iedge : iedges) {
-          const auto source = m_data.source(iedge);
-          const auto target = m_data.target(iedge);
-          const auto s = m_data.to_2d(support_plane_idx, source);
-          const auto t = m_data.to_2d(support_plane_idx, target);
-          const Segment_2 segment(s, t); Point_2 inter;
-          if (KSR::intersection(query_segment, segment, inter)) {
-            found_intersection = true;
-            break;
-          }
-        }
-        if (!found_intersection) {
-          to_be_merged.push_back(std::make_pair(pface_i, pface_j));
-        }
-      }
-    }
-    if (is_debug) std::cout << "pairs to be merged: " << to_be_merged.size() << std::endl;
-    CGAL_assertion(to_be_merged.size() == 1);
-
-    if (is_debug) {
-      std::cout << "merge: " <<
-      m_data.str(to_be_merged[0].first) << " + " <<
-      m_data.str(to_be_merged[0].second) << std::endl;
-      std::cout << "centroid i: " << m_data.centroid_of_pface(to_be_merged[0].first) << std::endl;
-      std::cout << "centroid j: " << m_data.centroid_of_pface(to_be_merged[0].second) << std::endl;
-    }
-
-    CGAL_assertion_msg(false,
-    "TODO: POLYGON SPLITTER, MERGE COPLANAR PFACES!");
-  }
-
-  void add_bisectors(
-    const std::vector< std::vector<Point_2> >& original_faces,
-    const KSR::size_t support_plane_idx) {
-
-    CGAL_assertion_msg(false,
-    "WARNING: WHEN ADDING BISECTORS, WE SHOULD ALSO CHANGE THE KINETIC CODE! NOT FINISHED!");
-
-    const bool is_debug = false;
-    CGAL_assertion(support_plane_idx >= 6);
-
-    if (is_debug) {
-      std::cout << std::endl << "support plane idx: " << support_plane_idx << std::endl;
-      std::cout << "dumping support plane ... ";
-      dump(true, 0, support_plane_idx, "support-plane-" +
-      std::to_string(support_plane_idx) + ".ply");
-      std::cout << "done" << std::endl;
-    }
-
-    std::vector<Weighted_point> wps;
-    create_weighted_points(original_faces, support_plane_idx, wps);
-
-    std::vector<RVertex_handle> vhs;
-    vhs.reserve(original_faces.size());
-    Regular_triangulation triangulation;
-    create_regular_triangulation(wps, vhs, triangulation);
-    CGAL_assertion(triangulation.is_valid());
-    CGAL_assertion(vhs.size() == original_faces.size());
-
-    if (is_debug) {
-      std::cout << "dumping regular triangulation / power diagram ... ";
-      dump(triangulation, support_plane_idx);
-      std::cout << "done" << std::endl;
-
-      for (std::size_t i = 0; i < vhs.size(); ++i) {
-        std::cout << "vhs " << std::to_string(i) << ": "
-        << m_data.to_3d(support_plane_idx, vhs[i]->point().point()) << std::endl;
-      }
-    }
-
-    std::vector< std::pair<Segment_2, bool> > bisectors;
-    find_bisectors(is_debug, support_plane_idx, vhs, triangulation, bisectors);
-    CGAL_assertion(bisectors.size() > 0);
-    convert_bisectors_to_iedges(is_debug, support_plane_idx, wps, bisectors);
-
-    // CGAL_assertion_msg(false,
-    // "TODO: POLYGON SPLITTER, ADD BISECTORS!");
-  }
-
-  void create_weighted_points(
-    const std::vector< std::vector<Point_2> >& original_faces,
-    const KSR::size_t support_plane_idx,
-    std::vector<Weighted_point>& wps) const {
-
-    wps.clear();
-    wps.reserve(original_faces.size() + 4);
-    add_weighted_bbox(support_plane_idx, wps);
-    CGAL_assertion(wps.size() == 4);
-    add_weighted_input(original_faces, support_plane_idx, wps);
-    CGAL_assertion(wps.size() == original_faces.size() + 4);
-  }
-
-  void add_weighted_bbox(
-    const KSR::size_t support_plane_idx,
-    std::vector<Weighted_point>& wps) const {
+    std::vector<Point_2>& bbox) const {
 
     CGAL_assertion(support_plane_idx >= 6);
     const auto& iedges = m_data.iedges(support_plane_idx);
@@ -863,176 +678,18 @@ private:
     }
     CGAL_assertion(points.size() == iedges.size() * 2);
 
-    const auto bbox = CGAL::bbox_2(points.begin(), points.end());
-    const Weighted_point wp1(Point_2(bbox.xmin(), bbox.ymin()), FT(0));
-    const Weighted_point wp2(Point_2(bbox.xmax(), bbox.ymin()), FT(0));
-    const Weighted_point wp3(Point_2(bbox.xmax(), bbox.ymax()), FT(0));
-    const Weighted_point wp4(Point_2(bbox.xmin(), bbox.ymax()), FT(0));
-    wps.push_back(wp1);
-    wps.push_back(wp2);
-    wps.push_back(wp3);
-    wps.push_back(wp4);
-  }
+    const auto box = CGAL::bbox_2(points.begin(), points.end());
+    const Point_2 p1(box.xmin(), box.ymin());
+    const Point_2 p2(box.xmax(), box.ymin());
+    const Point_2 p3(box.xmax(), box.ymax());
+    const Point_2 p4(box.xmin(), box.ymax());
 
-  void add_weighted_input(
-    const std::vector< std::vector<Point_2> >& original_faces,
-    const KSR::size_t support_plane_idx,
-    std::vector<Weighted_point>& wps) const {
-
-    for (const auto& original_face : original_faces) {
-      const auto centroid = CGAL::centroid(
-        original_face.begin(), original_face.end());
-
-      FT max_sq_dist = -FT(1);
-      for (const auto& p : original_face) {
-        const FT sq_dist = CGAL::squared_distance(p, centroid);
-        max_sq_dist = (CGAL::max)(sq_dist, max_sq_dist);
-      }
-      CGAL_assertion(max_sq_dist > FT(0));
-      const FT weight = static_cast<FT>(CGAL::sqrt(CGAL::to_double(max_sq_dist)));
-      const Weighted_point wp(centroid, weight);
-      wps.push_back(wp);
-    }
-  }
-
-  void create_regular_triangulation(
-    const std::vector<Weighted_point>& wps,
-    std::vector<RVertex_handle>& vhs,
-    Regular_triangulation& triangulation) const {
-
-    for (std::size_t i = 0; i < wps.size(); ++i) {
-      const auto& wp = wps[i];
-      if (i < 4) { triangulation.insert(wp); }
-      else {
-        const auto vh = triangulation.insert(wp);
-        vhs.push_back(vh);
-      }
-    }
-  }
-
-  void find_bisectors(
-    const bool is_debug,
-    const KSR::size_t support_plane_idx,
-    const std::vector<RVertex_handle>& vhs,
-    const Regular_triangulation& triangulation,
-    std::vector< std::pair<Segment_2, bool> >& bisectors) const {
-
-    if (vhs.size() > 2) {
-      CGAL_assertion_msg(false, "TODO: FIND MULTIPLE BISECTORS!");
-      // The way to find multiple bisectors is actually almost the same as for
-      // one bisector. We first find all pairs of unique vhs connections:
-      // like vhs[0] -> vhs[1], vhs[1] -> vhs[3], vhs[3] -> vhs[2], vhs[2] -> vhs[1], etc.
-      // then we extract the corresponding segments + possibly two rays, which
-      // may intersect the bbox boundaries. We transform these rays into segments.
-      // That is we are going to have several interior bisectors and two boundary bisectors.
-    }
-
-    bisectors.clear();
-    const auto edge = find_bisecting_edge(
-      is_debug, vhs[0], vhs[1], support_plane_idx, vhs, triangulation);
-    const auto object = triangulation.dual(edge);
-    Segment_2 bisector;
-    if (CGAL::assign(bisector, object)) {
-      if (is_debug) {
-        std::cout << "found bisector: 2 " <<
-        m_data.to_3d(support_plane_idx, bisector.source()) << " " <<
-        m_data.to_3d(support_plane_idx, bisector.target()) << std::endl;
-      }
-    } else {
-      CGAL_assertion_msg(false, "TODO: ADD OTHER TYPES: RAY/LINE!");
-    }
-    bisectors.push_back(std::make_pair(bisector, true));
-  }
-
-  const REdge find_bisecting_edge(
-    const bool is_debug,
-    const RVertex_handle& vh0,
-    const RVertex_handle& vh1,
-    const KSR::size_t support_plane_idx,
-    const std::vector<RVertex_handle>& vhs,
-    const Regular_triangulation& triangulation) const {
-
-    auto curr = triangulation.incident_edges(vh0);
-    const auto end = curr;
-    do {
-      const auto fh = curr->first;
-      const auto id = curr->second;
-      const auto ip = (id + 1) % 3;
-      const auto im = (id + 2) % 3;
-      const auto& p = fh->vertex(ip)->point().point();
-      const auto& q = fh->vertex(im)->point().point();
-      // std::cout << "ip, p: " << m_data.to_3d(support_plane_idx, p) << std::endl;
-      // std::cout << "im, q: " << m_data.to_3d(support_plane_idx, q) << std::endl;
-
-      CGAL_assertion(
-        fh->vertex(ip) == vh0 || fh->vertex(im) == vh0);
-      if (fh->vertex(ip) == vh1) {
-        if (is_debug) {
-          std::cout << "ip, found connecting edge: 2 " <<
-            m_data.to_3d(support_plane_idx, q) << " " <<
-            m_data.to_3d(support_plane_idx, p) << std::endl;
-        }
-        return *curr;
-      }
-
-      if (fh->vertex(im) == vh1) {
-        if (is_debug) {
-          std::cout << "im, found connecting edge: 2 " <<
-            m_data.to_3d(support_plane_idx, p) << " " <<
-            m_data.to_3d(support_plane_idx, q) << std::endl;
-        }
-        return *curr;
-      }
-
-      ++curr;
-    } while (curr != end);
-    CGAL_assertion_msg(curr == end, "ERROR: THE BISECTING EDGE IS NOT FOUND!");
-    return *curr;
-  }
-
-  void convert_bisectors_to_iedges(
-    const bool is_debug,
-    const KSR::size_t support_plane_idx,
-    const std::vector<Weighted_point>& wps,
-    const std::vector< std::pair<Segment_2, bool> >& bisectors) {
-
-    CGAL_assertion(bisectors.size() > 0);
-    if (bisectors.size() > 1) {
-      CGAL_assertion_msg(false, "TODO: HANDLE MULTIPLE BISECTORS!");
-      // Should be more or less the same code as for one bisector!
-    }
-
-    const bool is_boundary_bisector = bisectors[0].second;
-    CGAL_assertion(is_boundary_bisector);
-    const auto& bisector = bisectors[0].first;
-    std::vector<Point_2> inter_points;
-    for (std::size_t i = 0; i < 4; ++i) {
-      const std::size_t ip = (i + 1) % 4;
-      const Segment_2 bbox_edge(wps[i].point(), wps[ip].point());
-      Point_2 inter_point;
-      if (KSR::intersection(bisector, bbox_edge, inter_point)) {
-        inter_points.push_back(inter_point);
-      }
-    }
-    CGAL_assertion(inter_points.size() == 2);
-    if (is_debug) {
-      std::cout << "found iedge: 2 " <<
-        m_data.to_3d(support_plane_idx, inter_points[0]) << " " <<
-        m_data.to_3d(support_plane_idx, inter_points[1]) << std::endl;
-    }
-
-    const auto iv0 = m_data.igraph().
-      add_vertex(m_data.to_3d(support_plane_idx, inter_points[0])).first;
-    const auto iv1 = m_data.igraph().
-      add_vertex(m_data.to_3d(support_plane_idx, inter_points[1])).first;
-
-    const auto pair = m_data.igraph().add_edge(iv0, iv1, support_plane_idx);
-    const auto& iedge = pair.first;
-    const bool is_inserted = pair.second;
-    if (is_inserted) {
-      m_data.igraph().set_line(iedge, m_data.igraph().add_line());
-    }
-    m_data.support_plane(support_plane_idx).iedges().insert(iedge);
+    bbox.clear();
+    bbox.reserve(4);
+    bbox.push_back(p1);
+    bbox.push_back(p2);
+    bbox.push_back(p3);
+    bbox.push_back(p4);
   }
 
   void dump(
@@ -1119,56 +776,6 @@ private:
     KSR_3::Saver<Kernel> saver;
     saver.export_polygon_soup_3(polygons, file_name);
   }
-
-  void dump(
-    const Regular_triangulation& triangulation,
-    const KSR::size_t support_plane_idx) {
-
-    Mesh_3 mesh;
-    std::map<RVertex_handle, Vertex_index> map_v2i;
-    for (auto vit = triangulation.finite_vertices_begin();
-    vit != triangulation.finite_vertices_end(); ++vit) {
-      const auto wp = vit->point();
-      const Point_2 point(wp.x(), wp.y());
-      map_v2i.insert(std::make_pair(
-        vit, mesh.add_vertex(m_data.support_plane(support_plane_idx).to_3d(point))));
-    }
-
-    for (auto fit = triangulation.finite_faces_begin();
-    fit != triangulation.finite_faces_end(); ++fit) {
-      std::array<Vertex_index, 3> vertices;
-      for (std::size_t i = 0; i < 3; ++i) {
-        vertices[i] = map_v2i[fit->vertex(i)];
-      }
-      mesh.add_face(vertices);
-    }
-
-    std::ofstream out_cdt("power-cdt.ply");
-    out_cdt.precision(20);
-    CGAL::write_ply(out_cdt, mesh);
-    out_cdt.close();
-
-    std::ofstream out_diagram("power-diagram.polylines.txt");
-    out_diagram.precision(20);
-
-    for(auto eit = triangulation.finite_edges_begin();
-    eit != triangulation.finite_edges_end(); ++eit) {
-      const auto object = triangulation.dual(eit);
-
-      // typename Kernel::Ray_2  ray;
-      // typename Kernel::Line_2 line;
-
-      typename Kernel::Segment_2 segment;
-      if (CGAL::assign(segment, object)) {
-        out_diagram << "2 " <<
-        m_data.to_3d(support_plane_idx, segment.source()) << " " <<
-        m_data.to_3d(support_plane_idx, segment.target()) << std::endl;
-      }
-      // if (CGAL::assign(ray, object)) out_diagram << "2 " << r << std::endl;
-      // if (CGAL::assign(line, object)) out_diagram << "2 " << l << std::endl;
-    }
-    out_diagram.close();
-  }
 };
 
 } // namespace KSR_3
diff --git a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Support_plane.h b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Support_plane.h
index ca2e034c11d..a2e60fc1e5c 100644
--- a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Support_plane.h
+++ b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Support_plane.h
@@ -122,7 +122,7 @@ public:
       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: POLYGON HAS FLAT BBOX!");
+    CGAL_assertion_msg(normal != CGAL::NULL_VECTOR, "ERROR: BBOX IS FLAT!");
 
     m_data->k = 0;
     m_data->plane = Plane_3(points[0], KSR::normalize(normal));
diff --git a/Kinetic_shape_reconstruction/include/CGAL/Kinetic_shape_reconstruction_3.h b/Kinetic_shape_reconstruction/include/CGAL/Kinetic_shape_reconstruction_3.h
index 022d2e355f8..67d79af596e 100644
--- a/Kinetic_shape_reconstruction/include/CGAL/Kinetic_shape_reconstruction_3.h
+++ b/Kinetic_shape_reconstruction/include/CGAL/Kinetic_shape_reconstruction_3.h
@@ -82,7 +82,7 @@ private:
 public:
   Kinetic_shape_reconstruction_3(
     const bool verbose = true,
-    const bool debug = false) :
+    const bool debug   = false) :
   m_debug(debug),
   m_verbose(verbose),
   m_queue(m_debug),
@@ -92,51 +92,50 @@ public:
   m_data(m_debug)
   { }
 
-  // TODO: Use named parameters here!
+  // TODO: USE NAMED PARAMETERS!
   template<
   typename InputRange,
   typename PolygonMap>
   const bool partition(
     const InputRange& input_range,
     const PolygonMap polygon_map,
-    const unsigned int k = 1,
-    const unsigned int n = 0,
-    const double enlarge_bbox_ratio = 1.1,
-    const bool reorient = false) {
+    unsigned int k = 1,
+    unsigned int n = 0,
+    double enlarge_bbox_ratio = 1.1,
+    bool reorient = false) {
 
     std::cout.precision(20);
     if (input_range.size() == 0) {
-      CGAL_warning_msg(input_range.size() != 0,
-      "WARNING: YOUR INPUT IS EMPTY. RETURN WITH NO CHANGE!");
-      return false;
-    }
-
-    if (k == 0) {
-      CGAL_warning_msg(k != 0,
-      "WARNING: YOU SET K TO 0. THE VALID VALUES ARE {1,2,...}. RETURN WITH NO CHANGE!");
+      CGAL_warning_msg(input_range.size() > 0,
+      "WARNING: YOUR INPUT IS EMPTY! RETURN WITH NO CHANGE!");
       return false;
     }
 
     if (n != 0) {
-      CGAL_assertion_msg(n == 0,
-      "TODO: IMPLEMENT KINETIC SUBDIVISION!");
-
+      CGAL_assertion_msg(false, "TODO: IMPLEMENT KINETIC SUBDIVISION!");
       if (n > 3) {
         CGAL_warning_msg(n <= 3,
-        "WARNING: DOES IT MAKE SENSE TO HAVE MORE THAN 64 INPUT BLOCKS? RETURN WITH NO CHANGE!");
-        return false;
+        "WARNING: DOES IT MAKE SENSE TO HAVE MORE THAN 64 INPUT BLOCKS? SETTING N TO 3!");
+        n = 3;
       }
     }
 
     if (enlarge_bbox_ratio < 1.0) {
       CGAL_warning_msg(enlarge_bbox_ratio >= 1.0,
-      "WARNING: YOU SET ENLARGE_BBOX_RATIO < 1.0. THE VALID RANGE IS [1.0, +INF). RETURN WITH NO CHANGE!");
-      return false;
+      "WARNING: YOU SET ENLARGE_BBOX_RATIO < 1.0! THE VALID RANGE IS [1.0, +INF). SETTING TO 1.0!");
+      enlarge_bbox_ratio = 1.0;
     }
 
     const FT time_step = static_cast<FT>(m_initializer.initialize(
       input_range, polygon_map, k, enlarge_bbox_ratio, reorient));
     m_initializer.convert(m_data);
+    m_data.check_integrity();
+
+    if (k == 0) {
+      CGAL_warning_msg(k > 0,
+      "WARNING: YOU SET K TO 0! THAT MEANS NO PROPAGATION! THE VALID VALUES ARE {1,2,...}. INTERSECT AND RETURN!");
+      return false;
+    }
 
     // if (m_verbose) {
     //   std::cout << std::endl << "POLYGON SPLITTER SUCCESS!" << std::endl << std::endl;
@@ -803,7 +802,8 @@ private:
 
     const bool is_event_found = false;
     return is_event_found;
-    CGAL_assertion_msg(false, "TODO: ADD THIS EXTRA TYPE OF EVENT!");
+    CGAL_assertion_msg(false,
+    "TODO: ADD PVERTEX TO IVERTEX UNCONSTRAINED EVENT!");
   }
 
   const std::size_t run(
@@ -904,7 +904,7 @@ private:
     const Event&   /* event   */) {
 
     CGAL_assertion_msg(false,
-    "TODO: ADD CASE TWO CONSTRAINED PVERTICES MEET!");
+    "TODO: IMPLEMENT TWO CONSTRAINED PVERTICES MEET EVENT!");
   }
 
   void apply_event_two_unconstrained_pvertices_meet(
@@ -1192,7 +1192,7 @@ private:
     const Event&   /* event   */) {
 
     CGAL_assertion_msg(false,
-    "TODO: ADD UNCONSTRAINED PVERTEX MEETS IVERTEX EVENT!");
+    "TODO: IMPLEMENT UNCONSTRAINED PVERTEX MEETS IVERTEX EVENT!");
   }
 
   void remove_events(const IEdge& iedge, const KSR::size_t support_plane_idx) {
diff --git a/Kinetic_shape_reconstruction/test/Kinetic_shape_reconstruction/CMakeLists.txt b/Kinetic_shape_reconstruction/test/Kinetic_shape_reconstruction/CMakeLists.txt
index 3a11f2f1a33..433fc4ad82b 100644
--- a/Kinetic_shape_reconstruction/test/Kinetic_shape_reconstruction/CMakeLists.txt
+++ b/Kinetic_shape_reconstruction/test/Kinetic_shape_reconstruction/CMakeLists.txt
@@ -19,21 +19,30 @@ if(CGAL_FOUND)
   if(Boost_FOUND)
     message(STATUS "Found Boost")
 
-    set(targets
-    # kinetic_2d_stress_test
-    kinetic_3d_test_all
-    )
+    find_package(Eigen3 3.1.0 REQUIRED)
+    if(Eigen3_FOUND)
+      message(STATUS "Found Eigen")
 
-    set(project_linked_libraries)
-    set(project_compilation_definitions)
+      include(CGAL_Eigen_support)
 
-    foreach(target ${targets})
-      create_single_source_cgal_program("${target}.cpp")
-      if(TARGET ${target})
-        target_link_libraries(${target} PUBLIC ${project_linked_libraries})
-        target_compile_definitions(${target} PUBLIC ${project_compilation_definitions})
-      endif()
-    endforeach()
+      set(targets
+      # kinetic_2d_stress_test
+      kinetic_3d_test_all
+      )
+
+      set(project_linked_libraries)
+      set(project_compilation_definitions)
+
+      foreach(target ${targets})
+        create_single_source_cgal_program("${target}.cpp")
+        if(TARGET ${target})
+          target_link_libraries(${target} PUBLIC ${project_linked_libraries} CGAL::Eigen_support)
+          target_compile_definitions(${target} PUBLIC ${project_compilation_definitions})
+        endif()
+      endforeach()
+    else()
+      message(ERROR "This program requires the Eigen library, and will not be compiled.")
+    endif()
   else()
     message(ERROR "This program requires the Boost library, and will not be compiled.")
   endif()
diff --git a/Kinetic_shape_reconstruction/test/Kinetic_shape_reconstruction/kinetic_3d_test_all.cpp b/Kinetic_shape_reconstruction/test/Kinetic_shape_reconstruction/kinetic_3d_test_all.cpp
index 1e3087190c9..10b4dd6a722 100644
--- a/Kinetic_shape_reconstruction/test/Kinetic_shape_reconstruction/kinetic_3d_test_all.cpp
+++ b/Kinetic_shape_reconstruction/test/Kinetic_shape_reconstruction/kinetic_3d_test_all.cpp
@@ -139,11 +139,8 @@ int main (const int argc, const char** argv) {
   assert(run_test("data/edge-case-test/test-2-polygons.off" , ks, num_iters, num_tests)); // edge touch
   assert(run_test("data/edge-case-test/test-4-polygons.off" , ks, num_iters, num_tests)); // edge touch / 2 coplanar
   assert(run_test("data/edge-case-test/test-5-polygons.off" , ks, num_iters, num_tests)); // edge touch / vertex touch / 2 coplanar
-
-  // std::vector<unsigned int> ts;
-  // ts.push_back(1); ts.push_back(2); ts.push_back(4);
-  // ts.push_back(5); ts.push_back(6); ts.push_back(100);
   assert(run_test("data/edge-case-test/test-20-polygons.off", ks, num_iters, num_tests)); // 2 overlap and coplanar
+  assert(run_test("data/edge-case-test/test-flat-bbox.off"  , ks, num_iters, num_tests)); // flat bbox / 2 coplanar
 
   std::cout << std::endl << "--OUTPUT STATS:" << std::endl;
   std::cout << "* number of iterations per test: " << num_iters << std::endl;