diff --git a/Snap_rounding_2/examples/Snap_rounding_2/float_snap_rounding.cpp b/Snap_rounding_2/examples/Snap_rounding_2/float_snap_rounding.cpp
index 4c4aecd6b3a..4866de1edf0 100644
--- a/Snap_rounding_2/examples/Snap_rounding_2/float_snap_rounding.cpp
+++ b/Snap_rounding_2/examples/Snap_rounding_2/float_snap_rounding.cpp
@@ -9,13 +9,14 @@ typedef CGAL::Arr_segment_traits_2<Kernel>                      Traits_2;
 typedef Traits_2::Curve_2                                       Segment_2;
 typedef Kernel::Point_2                          Point_2;
 typedef Kernel::FT                               FT;
+typedef std::vector<Point_2 > Polyline_2;
 
 int main()
 {
   std::vector<Point_2> pts;
   std::vector< Segment_2 > segs;
 
-  Kernel::FT e(std::pow(2, -60));
+  FT e(std::pow(2, -60));
 
   segs.emplace_back(Point_2(1-e, 1), Point_2(-1-e, -1+2*e));
   segs.emplace_back(Point_2(e/2, e/2), Point_2(1, -1));
@@ -31,7 +32,7 @@ int main()
   }
   std::cout << "\n\n" << std::endl;
 
-  std::vector< std::vector<Point_2 > > out;
+  std::vector< Polyline_2 > out;
   CGAL::double_snap_rounding_2(segs.begin(), segs.end(), out);
 
   std::cout << "Output" << std::endl;
diff --git a/Snap_rounding_2/include/CGAL/Float_snap_rounding_2.h b/Snap_rounding_2/include/CGAL/Float_snap_rounding_2.h
index 4f60d5a9b13..7720e709dd9 100644
--- a/Snap_rounding_2/include/CGAL/Float_snap_rounding_2.h
+++ b/Snap_rounding_2/include/CGAL/Float_snap_rounding_2.h
@@ -125,10 +125,12 @@ void double_snap_rounding_2_disjoint(PointsRange &pts,
     return res==SMALLER;
   };
 
-
+#ifdef DOUBLE_2D_SNAP_VERBOSE
+  std::cout << "Sort the input points" << std::endl;
+#endif
   // Compute the order of the points along the 2 axis
   // Sorted the points may perform exact computations and thus refine the intervals of the coordinates values
-  // This refine ensures that the order of the points will be preserved by the rounding
+  // This refine ensures that the order of the points will be preserved when rounded
   using Iterator_set_x = typename std::set<size_t, decltype(comp_by_x_first)>::iterator;
   using Iterator_set_y = typename std::set<size_t, decltype(comp_by_y_first)>::iterator;
   std::set<size_t, decltype(comp_by_x_first)> p_sort_by_x(comp_by_x_first);
@@ -138,11 +140,6 @@ void double_snap_rounding_2_disjoint(PointsRange &pts,
     p_sort_by_x.insert(i);
     p_sort_by_y.insert(i);
   }
-  // std::vector<size_t> p_sort_by_y(pts.size(),0);
-  // std::iota(p_sort_by_x.begin(),p_sort_by_x.end(),0);
-  // std::iota(p_sort_by_y.begin(),p_sort_by_y.end(),0);
-  // std::sort(p_sort_by_x.begin(),p_sort_by_x.end(),comp_by_x);
-  // std::sort(p_sort_by_y.begin(),p_sort_by_y.end(),comp_by_y);
 
   //Kd-Tree to exhibits pairs
   std::vector<PBox> points_boxes;
@@ -152,8 +149,10 @@ void double_snap_rounding_2_disjoint(PointsRange &pts,
   for(size_t i=0; i<polylines.size(); ++i)
     segs_boxes.emplace_back(pts[polylines[i][0]].bbox()+pts[polylines[i][1]].bbox(),i);
 
+  // Bound the maximum squared distance between a point and its rounded value
   auto round_bound=[](Point_2& p){
-    return std::pow(p.x().approx().sup()-p.x().approx().inf(),2)+std::pow(p.y().approx().sup()-p.y().approx().inf(), 2);
+    return std::pow(std::nextafter(p.x().approx().sup(),std::numeric_limits<double>::infinity())-p.x().approx().inf(),2)+
+           std::pow(std::nextafter(p.y().approx().sup(),std::numeric_limits<double>::infinity())-p.y().approx().inf(),2);
   };
   auto callback=[&](PBox &bp, SBox &bseg){
     size_t pi=bp.index();
@@ -161,35 +160,54 @@ void double_snap_rounding_2_disjoint(PointsRange &pts,
     size_t si1=polylines[bseg.index()][0];
     size_t si2=polylines[bseg.index()][1];
 
-    // the point is a vertex of the segment
     if((pi==si1) || (pi==si2))
       return;
 
     Point_2& p= pts[bp.index()];
     Segment_2 seg(pts[polylines[bseg.index()][0]], pts[polylines[bseg.index()][1]]);
 
-    double round_bound_s=(std::max)(round_bound(pts[si1]), round_bound(pts[si2]));
+    // (A+B)^2 <= 4*max(A^2,B^2) and we take some margin
+    double bound=16*(std::max)(round_bound(pts[pi]), (std::max)(round_bound(pts[si1]), round_bound(pts[si2])));
 
-    if(possibly(Kernel().compare_squared_distance_2_object()(p, seg, round_bound_s)!=CGAL::LARGER)){
+    // If the segment is closed to the vertex and we subdivide it at same x coordinate of that vertex
+    if(possibly(Kernel().compare_squared_distance_2_object()(p, seg, bound)!=CGAL::LARGER) &&
+       compare(seg.source().x(),p.x())!=compare(seg.target().x(),p.x()))
+    {
       pts.emplace_back(p.x(), seg.supporting_line().y_at_x(p.x()));
       auto pair=p_sort_by_x.insert(pts.size()-1);
       if(pair.second)
         p_sort_by_y.insert(pts.size()-1);
       else
-        pts.pop_back();
+        pts.pop_back(); // Remove the new point if it is already exist
       polylines[si].emplace_back(*pair.first);
     }
   };
 
-  CGAL::box_intersection_d(points_boxes.begin(), points_boxes.end(), segs_boxes.begin(), segs_boxes.end(), callback);
+#ifdef DOUBLE_2D_SNAP_VERBOSE
+  std::cout << "Exhibit pairs of possible intersections" << std::endl;
+#endif
 
-  //sort new vertices
+  do{
+    size_t size_before=pts.size();
+    CGAL::box_intersection_d(points_boxes.begin(), points_boxes.end(), segs_boxes.begin(), segs_boxes.end(), callback);
+    points_boxes.clear();
+    // The new vertices may intersect a segment when rounded, we repeat until they are not new vertices
+#ifdef DOUBLE_2D_SNAP_VERBOSE
+    std::cout << points_boxes.size()-size_before << " subdivisions performed" << std::endl;
+#endif
+    for(size_t i=size_before; i<pts.size(); ++i)
+      points_boxes.emplace_back(pts[i].bbox(),i);
+  } while(points_boxes.size()!=0);
+
+#ifdef DOUBLE_2D_SNAP_VERBOSE
+  std::cout << "Form the polylines" << std::endl;
+#endif
   for(auto &polyline: polylines)
   {
     if(polyline.size()==2)
       continue;
 
-    //Sort the points on the polyline along the original vector
+    // Sort the subdivision points on the polyline along the original vector
     Vector_2 ref(pts[polyline[0]], pts[polyline[1]]);
     auto sort_along_ref=[&](size_t pi, size_t qi){
       Vector_2 v(pts[pi], pts[qi]);
@@ -197,10 +215,15 @@ void double_snap_rounding_2_disjoint(PointsRange &pts,
         return is_positive(v.y()*ref.y());
       return is_positive(v.x()*ref.x());
     };
+    size_t ps=polyline[0];
+    size_t pt=polyline[1];
     std::sort(polyline.begin(), polyline.end(), sort_along_ref);
+    CGAL_assertion((polyline[0]==ps) && (polyline[polyline.size()-1]==pt));
   }
 
-  //round
+#ifdef DOUBLE_2D_SNAP_VERBOSE
+  std::cout << "Round" << std::endl;
+#endif
   for(auto &p: pts)
     p=Point_2(to_double(p.x()), to_double(p.y()));
 
@@ -208,7 +231,9 @@ void double_snap_rounding_2_disjoint(PointsRange &pts,
   CGAL_assertion(std::is_sorted(p_sort_by_x.begin(),p_sort_by_x.end(),comp_by_x));
   CGAL_assertion(std::is_sorted(p_sort_by_y.begin(),p_sort_by_y.end(),comp_by_y));
 
-  //remove duplicate_points
+#ifdef DOUBLE_2D_SNAP_VERBOSE
+  std::cout << "Remove duplicate points" << std::endl;
+#endif
   std::vector< size_t > unique_points(p_sort_by_x.begin(),p_sort_by_x.end());
   std::sort(unique_points.begin(),unique_points.end(),comp_by_x_first);
   std::vector<Point_2> new_pts;
@@ -220,7 +245,6 @@ void double_snap_rounding_2_disjoint(PointsRange &pts,
   }
 
   std::swap(pts, new_pts);
-  // Update the polylines by remapping the old indices to new indices
   for (auto& polyline : polylines) {
     std::vector<size_t> updated_polyline;
     for (size_t i=0; i<polyline.size(); ++i) {
@@ -232,7 +256,11 @@ void double_snap_rounding_2_disjoint(PointsRange &pts,
     std::swap(polyline, updated_polyline);
   }
 
-  //Vertices can be on vertical or horizontal segment, we repair this
+#ifdef DOUBLE_2D_SNAP_VERBOSE
+  // The algorithm prevents the a vertex that goes through a segment but a vertex may lie on an horizontal/vertical segments after rounding
+  std::cout << "Subdivide horizontal and vertical segments with vertices on them" << std::endl;
+#endif
+  //The order may have changed (Example: (1,1)<(1,2)<(1+e,1))
   p_sort_by_x.clear();
   p_sort_by_y.clear();
   for(size_t i=0; i!=pts.size(); ++i)
@@ -240,16 +268,14 @@ void double_snap_rounding_2_disjoint(PointsRange &pts,
     p_sort_by_x.insert(i);
     p_sort_by_y.insert(i);
   }
-  // std::vector< size_t > vec_sort_by_x_first(p_sort_by_x.begin(),p_sort_by_x.end());
-  // std::vector< size_t > vec_sort_by_y_first(p_sort_by_y.begin(),p_sort_by_y.end());
-  // std::sort(vec_sort_by_x_first.begin(),vec_sort_by_x_first.end(),comp_by_x_first);
-  // std::sort(vec_sort_by_y_first.begin(),vec_sort_by_y_first.end(),comp_by_y_first);
+
   for(auto &poly: polylines){
     std::vector<size_t> updated_polyline;
     updated_polyline.push_back(poly[0]);
     for(size_t i=1; i!=poly.size(); ++i){
       if(pts[poly[i-1]].x()==pts[poly[i]].x()){
         Iterator_set_x start, end;
+        // Get all vertices between the two endpoints along x order
         if(comp_by_x_first(poly[i-1],poly[i])){
           start=p_sort_by_x.upper_bound(poly[i-1]);
           end=p_sort_by_x.lower_bound(poly[i]);
@@ -257,12 +283,14 @@ void double_snap_rounding_2_disjoint(PointsRange &pts,
           start=p_sort_by_x.upper_bound(poly[i]);
           end=p_sort_by_x.lower_bound(poly[i-1]);
         }
+        // Add all endpoints between them to the polyline
         for(auto it=start; it!=end; ++it){
           updated_polyline.push_back(*it);
         }
       }
       if(pts[poly[i-1]].y()==pts[poly[i]].y()){
         Iterator_set_y start, end;
+        // Get all vertices between the two endpoints along y order
         if(comp_by_y_first(poly[i-1],poly[i])){
           start=p_sort_by_y.upper_bound(poly[i-1]);
           end=p_sort_by_y.lower_bound(poly[i]);
@@ -270,6 +298,7 @@ void double_snap_rounding_2_disjoint(PointsRange &pts,
           start=p_sort_by_y.upper_bound(poly[i]);
           end=p_sort_by_y.lower_bound(poly[i-1]);
         }
+        // Add all endpoints between them to the polyline
         for(auto it=start; it!=end; ++it){
           updated_polyline.push_back(*it);
         }
@@ -278,9 +307,12 @@ void double_snap_rounding_2_disjoint(PointsRange &pts,
     }
     std::swap(poly, updated_polyline);
   }
-  // compute_subcurves(input_begin, input_end, polylines);
 }
 
+
+/*
+TODO doc
+*/
 template <class InputIterator , class OutputContainer>
 typename OutputContainer::iterator double_snap_rounding_2(InputIterator  	input_begin,
 		                                                      InputIterator  	input_end,
@@ -290,9 +322,15 @@ typename OutputContainer::iterator double_snap_rounding_2(InputIterator  	input_
   using Polyline = std::remove_cv_t<typename std::iterator_traits<typename OutputContainer::iterator>::value_type>;
   using Point_2 = typename Default_arr_traits<Segment_2>::Traits::Point_2;
 
-  std::vector<Segment_2> segs;
+  std::vector<Segment_2> segs(input_begin, input_end);
+#ifdef DOUBLE_2D_SNAP_VERBOSE
+  std::cout << "Solved intersections" << std::endl;
+#endif
   compute_subcurves(input_begin, input_end, std::back_inserter(segs));
 
+#ifdef DOUBLE_2D_SNAP_VERBOSE
+  std::cout << "Change format to range of points and indexes" << std::endl;
+#endif
   std::set<Point_2> unique_point_set;
   std::map<Point_2, int> point_to_index;
   std::vector<Point_2> pts;
@@ -317,25 +355,103 @@ typename OutputContainer::iterator double_snap_rounding_2(InputIterator  	input_
     }
   }
 
-  for(InputIterator it=input_begin; it!=input_end; ++it)
+  for(InputIterator it=segs.begin(); it!=segs.end(); ++it)
   {
     size_t index1 = point_to_index[it->source()];
     size_t index2 = point_to_index[it->target()];
     polylines.push_back({index1, index2});
   }
 
+  // Main algorithm
   double_snap_rounding_2_disjoint(pts, polylines);
+
+#ifdef DOUBLE_2D_SNAP_VERBOSE
+  std::cout << "Build output" << std::endl;
+#endif
+
+  // Output polylines
+  output.clear();
   for(auto &poly: polylines){
     Polyline new_line;
-    for(size_t pi: poly){
-        new_line.push_back(pts[pi]);
-    }
+    for(size_t pi: poly)
+      new_line.push_back(pts[pi]);
     output.push_back(new_line);
   }
 
   return output.begin();
 }
 
+/*
+TODO doc
+*/
+template <class InputIterator , class OutputContainer>
+typename OutputContainer::iterator compute_snap_subcurves_2(InputIterator  	input_begin,
+		                                                  InputIterator  	input_end,
+		                                                  OutputContainer&  output)
+{
+  using Segment_2 = std::remove_cv_t<typename std::iterator_traits<InputIterator>::value_type>;
+  using Point_2 = typename Default_arr_traits<Segment_2>::Traits::Point_2;
+
+  std::vector<Segment_2> segs;
+#ifdef DOUBLE_2D_SNAP_VERBOSE
+  std::cout << "Solved intersections" << std::endl;
+#endif
+
+  compute_subcurves(input_begin, input_end, std::back_inserter(segs));
+
+#ifdef DOUBLE_2D_SNAP_VERBOSE
+  std::cout << "Change format to range of points and indexes" << std::endl;
+#endif
+  std::set<Point_2> unique_point_set;
+  std::map<Point_2, int> point_to_index;
+  std::vector<Point_2> pts;
+  std::vector< std::vector< size_t> > polylines;
+
+  for(typename std::vector<Segment_2>::iterator it=segs.begin(); it!=segs.end(); ++it)
+  {
+    const Point_2& p1 = it->source();
+    const Point_2& p2 = it->target();
+
+    // Check and insert the endpoints if they are not already added
+    if (unique_point_set.find(p1) == unique_point_set.end()) {
+      unique_point_set.insert(p1);
+      pts.push_back(p1);
+      point_to_index[p1] = pts.size() - 1;
+    }
+    if (unique_point_set.find(p2) == unique_point_set.end()) {
+      unique_point_set.insert(p2);
+      pts.push_back(p2);
+      point_to_index[p2] = pts.size() - 1;
+    }
+  }
+
+  for(InputIterator it=segs.begin(); it!=segs.end(); ++it)
+  {
+    size_t index1 = point_to_index[it->source()];
+    size_t index2 = point_to_index[it->target()];
+    polylines.push_back({index1, index2});
+  }
+
+  // Main algorithm
+  double_snap_rounding_2_disjoint(pts, polylines);
+
+#ifdef DOUBLE_2D_SNAP_VERBOSE
+  std::cout << "Build output" << std::endl;
+#endif
+
+  // Output a range of segments while removing duplicate ones
+  std::set< std::pair<size_t,size_t> > set_out_segs;
+  output.clear();
+  for(auto &poly: polylines){
+    for(size_t i=1; i<poly.size(); ++i)
+      set_out_segs.emplace((std::min)(poly[i-1],poly[i]),(std::max)(poly[i-1],poly[i]));
+  }
+  for(auto &pair: set_out_segs)
+    output.emplace_back(pts[pair.first], pts[pair.second]);
+
+  return output.begin();
+}
+
 } //namespace CGAL
 
 #endif
diff --git a/Snap_rounding_2/test/Snap_rounding_2/CMakeLists.txt b/Snap_rounding_2/test/Snap_rounding_2/CMakeLists.txt
index 9add288c3ac..4a86b77fee0 100644
--- a/Snap_rounding_2/test/Snap_rounding_2/CMakeLists.txt
+++ b/Snap_rounding_2/test/Snap_rounding_2/CMakeLists.txt
@@ -6,7 +6,22 @@ project(Snap_rounding_2_Tests)
 
 find_package(CGAL REQUIRED)
 
-create_single_source_cgal_program(test_snap_rounding_2.cpp NO_TESTING)
+if(MSVC)
+  # Turn off a VC++ warning on a potential division by zero
+  # in Cartesian_kernel/include/CGAL/constructions/kernel_ftC3.h
+  # where CGAL_assume() does not help
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /wd4723")
+endif()
+
+# create a target per cppfile
+file(
+  GLOB cppfiles
+  RELATIVE ${CMAKE_CURRENT_SOURCE_DIR}
+  ${CMAKE_CURRENT_SOURCE_DIR}/*.cpp)
+foreach(cppfile ${cppfiles})
+  create_single_source_cgal_program("${cppfile}")
+endforeach()
+
 
 function(add_Snap_rounding_tests name)
   set(data_dir "data")