diff --git a/Arrangement_on_surface_2/include/CGAL/Arr_segment_traits_2.h b/Arrangement_on_surface_2/include/CGAL/Arr_segment_traits_2.h
index 13c7e27ba19..676b3c2bc1c 100644
--- a/Arrangement_on_surface_2/include/CGAL/Arr_segment_traits_2.h
+++ b/Arrangement_on_surface_2/include/CGAL/Arr_segment_traits_2.h
@@ -712,6 +712,11 @@ public:
       // Check if we have a single intersection point.
       const Point_2* ip = std::get_if<Point_2>(&*res);
       if (ip != nullptr) {
+        if(!(cv1.is_vertical() ?
+                       m_traits.is_in_y_range_2_object()(cv1, *ip) :
+                       m_traits.is_in_x_range_2_object()(cv1, *ip))){
+                        std::cout << cv1 << std::endl;
+                       }
         CGAL_assertion(cv1.is_vertical() ?
                        m_traits.is_in_y_range_2_object()(cv1, *ip) :
                        m_traits.is_in_x_range_2_object()(cv1, *ip));
diff --git a/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h b/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h
index 73c54858bf2..d99ac998b0e 100644
--- a/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h
+++ b/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h
@@ -406,6 +406,9 @@ CGAL_add_named_parameter(do_not_modify_geometry_t, do_not_modify_geometry, do_no
 CGAL_add_named_parameter_with_compatibility_ref_only(angles_param_t, angles_param, angles)
 CGAL_add_named_parameter(maximum_height_t, maximum_height, maximum_height)
 
+//List of named parameters used in Snap_rounding_2
+CGAL_add_named_parameter(compute_intersection_t, compute_intersection, compute_intersection)
+
 // List of named parameters used in the package 'Constrained_triangulation_3'
 CGAL_add_named_parameter(plc_face_id_t, plc_face_id, plc_face_id)
 CGAL_add_named_parameter(with_plc_face_id_t, with_plc_face_id, with_plc_face_id)
diff --git a/Snap_rounding_2/doc/Snap_rounding_2/Concepts/FloatSnapRoundingTraits_2.h b/Snap_rounding_2/doc/Snap_rounding_2/Concepts/FloatSnapRoundingTraits_2.h
new file mode 100644
index 00000000000..e067d0e6d18
--- /dev/null
+++ b/Snap_rounding_2/doc/Snap_rounding_2/Concepts/FloatSnapRoundingTraits_2.h
@@ -0,0 +1,228 @@
+
+/*!
+\ingroup PkgSnapRounding2Concepts
+\cgalConcept
+
+The concept `FloatSnapRoundingTraits_2` lists the set of requirements that must be fulfilled by
+an instance of the `Traits` template-parameter of
+the free functions of TODO.
+The list includes the nested types of the geometric primitives used in this class and
+some function object types for the required predicates on those primitives.
+
+\cgalRefines{AosTraits_2}
+
+\cgalHasModelsBegin
+\cgalHasModels{CGAL::Float_snap_rounding_traits_2<Kernel>}
+\cgalHasModelsEnd
+*/
+
+class FloatSnapRoundingTraits_2 {
+public:
+
+/// \name Types
+/// @{
+
+/*!
+The number type. This type must fulfill the requirements on
+`FieldNumberType`
+*/
+typedef unspecified_type FT;
+
+/*!
+Models the concept `ArrTraits::Point_2`.
+*/
+typedef unspecified_type Point_2;
+
+/*!
+Models the concept `ArrTraits::XMonotoneCurve_2`.
+*/
+typedef unspecified_type Segment_2;
+
+/// @}
+
+/// \name Functor Types
+/// @{
+
+/*!
+Models the concept `Kernel::ConstructSource_2`.
+*/
+typedef unspecified_type Construct_source_2;
+
+/*!
+Models the concept `Kernel::ConstructTarget_2`.
+*/
+typedef unspecified_type Construct_target_2;
+
+/*!
+Models the concept `Kernel::ConstructSegment_2`.
+*/
+typedef unspecified_type Construct_segment_2;
+
+/*!
+Models the concept `Kernel::LessXY_2`.
+*/
+typedef unspecified_type Less_xy_2;
+
+/*!
+Models the concept `Kernel::LessY_2`.
+*/
+typedef unspecified_type Less_y_2;
+
+/*!
+Models the concept `Kernel::Equal_2`.
+*/
+typedef unspecified_type Equal_2;
+
+/*!
+Models the concept `FSRTraits_2::ConstructRoundPoint_2`.
+*/
+typedef unspecified_type Construct_round_point_2;
+
+/*!
+Models the concept `FSRTraits_2::SquaredRoundBound_2`.
+*/
+typedef unspecified_type Squared_round_bound_2;
+
+/*!
+Models the concept `FSRTraits_2::ConverterToExact`.
+*/
+typedef unspecified_type Converter_to_exact;
+
+/*!
+Models the concept `FSRTraits_2::ConverterFromExact`.
+*/
+typedef unspecified_type Converter_from_exact;
+
+/// @}
+
+/// \name Accessing Functor Objects
+/// @{
+
+/*!
+
+*/
+Construct_source_2 construct_source_2_object();
+
+/*!
+
+*/
+Construct_target_2 construct_target_2_object();
+
+/*!
+
+*/
+Construct_segment_2 construct_segment_2_object();
+
+/*!
+
+*/
+Less_xy_2 less_xy_2_object();
+
+/*!
+
+*/
+Less_y_2 less_y_2_object();
+
+/*!
+
+*/
+Construct_round_point_2 construct_round_point_2_object();
+
+/*!
+
+*/
+Squared_round_bound_2 squared_round_bound_2_object();
+
+/*!
+
+*/
+Converter_to_exact converter_to_exact_object();
+
+/*!
+
+*/
+Converter_from_exact converter_from_exact_object();
+
+
+/// @}
+
+}; /* end FloatSnapRoundingTraits_2 */
+
+
+namespace FSRTraits_2{
+
+/*!
+  \ingroup PkgSnapRounding2Concepts
+  \cgalConcept
+  \cgalHasModelsBegin
+  \cgalHasModelsBare{\link FloatSnapRoundingTraits_2::Construct_round_point_2 `Float_snap_rounding_traits_2::Construct_round_point_2` \endlink}
+  \cgalHasModelsEnd
+*/
+class ConstructRoundPoint_2
+{
+  public:
+
+  /*!
+  Given a point, construct its rounded version
+  */
+  Point_2 operator()(Point_2 p);
+};
+
+/*!
+  \ingroup PkgSnapRounding2Concepts
+  \cgalConcept
+  \cgalHasModelsBegin
+  \cgalHasModelsBare{\link FloatSnapRoundingTraits_2::Squared_round_bound_2 `Float_snap_rounding_traits_2::Squared_round_bound_2` \endlink}
+  \cgalHasModelsEnd
+*/
+class SquaredRoundBound_2
+{
+  public:
+
+  /*!
+  Given a point, compute an upper bound of the squared distance between its exact coordinates and its rounded coordinates
+  */
+  double operator()(Point_2 p);
+};
+
+/*!
+  \ingroup PkgSnapRounding2Concepts
+  \cgalConcept
+  \cgalHasModelsBegin
+  \cgalHasModelsBare{\link FloatSnapRoundingTraits_2::Converter_to_exact `Float_snap_rounding_traits_2::Convertex_to_exact` \endlink}
+  \cgalHasModelsEnd
+*/
+class ConverterToExact
+{
+  public:
+  typedef unspecified_type InputPoint_2;
+  typedef unspecified_type InputSegment_2;
+  /*!
+  Convert input points (i.e. segments) into points (i.e. segments) of the type of the traits
+  */
+  Point_2 operator()(InputPoint_2 p);
+  Segment_2 operator()(InputSegment_2 p);
+};
+
+/*!
+  \ingroup PkgSnapRounding2Concepts
+  \cgalConcept
+  \cgalHasModelsBegin
+  \cgalHasModelsBare{\link FloatSnapRoundingTraits_2::Converter_from_exact `Float_snap_rounding_traits_2::Convertex_to_exact` \endlink}
+  \cgalHasModelsEnd
+*/
+class ConverterFromExact
+{
+  public:
+  typedef unspecified_type OutputPoint_2;
+  typedef unspecified_type OutputSegment_2;
+  /*!
+  Convert points (i.e. segments) of the type of the traits into points (i.e. segments) of the output
+  */
+  OutputPoint_2 operator()(Point_2 p);
+  OutputSegment_2 operator()(Segment_2 p);
+};
+
+
+
+} /* end of namespace SRTraits_2 */
diff --git a/Snap_rounding_2/doc/Snap_rounding_2/Doxyfile.in b/Snap_rounding_2/doc/Snap_rounding_2/Doxyfile.in
index 7c99d29c3b1..0e85668587a 100644
--- a/Snap_rounding_2/doc/Snap_rounding_2/Doxyfile.in
+++ b/Snap_rounding_2/doc/Snap_rounding_2/Doxyfile.in
@@ -1,3 +1,11 @@
 @INCLUDE = ${CGAL_DOC_PACKAGE_DEFAULTS}
 
 PROJECT_NAME = "CGAL ${CGAL_DOC_VERSION} - 2D Snap Rounding"
+
+# custom options for this package
+HIDE_UNDOC_CLASSES         =  true
+WARN_IF_UNDOCUMENTED       =  false
+
+
+INPUT += ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Float_snap_rounding_2.h \
+         ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Float_snap_rounding_traits_2.h
diff --git a/Snap_rounding_2/examples/Snap_rounding_2/CMakeLists.txt b/Snap_rounding_2/examples/Snap_rounding_2/CMakeLists.txt
index e09fa21f5d5..315389de4a2 100644
--- a/Snap_rounding_2/examples/Snap_rounding_2/CMakeLists.txt
+++ b/Snap_rounding_2/examples/Snap_rounding_2/CMakeLists.txt
@@ -14,3 +14,11 @@ file(
 foreach(cppfile ${cppfiles})
   create_single_source_cgal_program("${cppfile}")
 endforeach()
+
+find_package(TBB QUIET)
+include(CGAL_TBB_support)
+if(TARGET CGAL::TBB_support)
+  target_link_libraries(float_snap_rounding PRIVATE CGAL::TBB_support)
+else()
+  message(STATUS "NOTICE: Intel TBB was not found. Sequential code will be used.")
+endif()
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
new file mode 100644
index 00000000000..9ac611ca936
--- /dev/null
+++ b/Snap_rounding_2/examples/Snap_rounding_2/float_snap_rounding.cpp
@@ -0,0 +1,68 @@
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Float_snap_rounding_2.h>
+#include <CGAL/Arr_segment_traits_2.h>
+#include <CGAL/Random.h>
+
+typedef CGAL::Exact_predicates_exact_constructions_kernel       Kernel;
+typedef CGAL::Exact_predicates_inexact_constructions_kernel     Epick;
+typedef Kernel::Segment_2                                       Segment_2;
+typedef Kernel::Point_2                                         Point_2;
+typedef Kernel::FT                                              FT;
+typedef std::vector<Point_2 >                                   Polyline_2;
+
+int main(int argc, char *argv[])
+{
+  std::vector< Segment_2 > segs;
+
+  if(argc>1){
+    std::cout << "Read segments in " << argv[1] << std::endl;
+    std::ifstream in(argv[1]);
+
+    int n;
+    in >> n;
+    segs.reserve(n);
+    for (int i=0; i<n; ++i)
+    {
+      double x1,x2,y1,y2;
+      in >> x1 >> y1 >> x2 >> y2;
+      if(Point_2(x1,y1)!=Point_2(x2,y2))
+        segs.emplace_back(Point_2(x1, y1), Point_2(x2, y2));
+    }
+  } else {
+    //Example with a non-trivial rounding
+    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));
+    segs.emplace_back(Point_2(0, 2-e/2), Point_2(2, 0));
+    segs.emplace_back(Point_2(0, 2-e/2), Point_2(-2+e, -4));
+    segs.emplace_back(Point_2(-2, 2), Point_2(2, 2));
+    segs.emplace_back(Point_2(7, 7), Point_2(7+e, 7+e));
+    segs.emplace_back(Point_2(5, 7-e), Point_2(9, 7-e));
+  }
+
+  std::cout << "Computes the intersections and snaps the segments" << std::endl;
+  std::vector< Segment_2> out;
+  CGAL::compute_snapped_subcurves_2(segs.begin(), segs.end(), std::back_inserter(out));
+  std::cout << "Does the output intersect: " << CGAL::do_curves_intersect(out.begin(), out.end()) << std::endl;
+  std::cout << "Size of the output: " << out.size() << std::endl;
+
+  std::string out_path=(argc>2)?argv[2]:"out.segs";
+  std::cout << "Write the outputs in " << out_path << std::endl;
+  std::ofstream outf(out_path);
+  CGAL::Random r;
+
+  outf << std::setprecision(17);
+  outf << out.size() << std::endl;
+  for (auto &seg: out)
+  {
+    double x1 = CGAL::to_double(seg.source().x());
+    double y1 = CGAL::to_double(seg.source().y());
+    double x2 = CGAL::to_double(seg.target().x());
+    double y2 = CGAL::to_double(seg.target().y());
+    outf << x1 << " " << y1 << " " << x2 << " " << y2 << std::endl;
+  }
+
+  return 0;
+}
diff --git a/Snap_rounding_2/include/CGAL/Float_snap_rounding_2.h b/Snap_rounding_2/include/CGAL/Float_snap_rounding_2.h
new file mode 100644
index 00000000000..66d1630fdcc
--- /dev/null
+++ b/Snap_rounding_2/include/CGAL/Float_snap_rounding_2.h
@@ -0,0 +1,830 @@
+// Copyright (c) 2025 Geometry Factory.
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+//
+// author(s)     : Léo Valque
+
+#ifndef CGAL_FLOAT_SNAP_ROUNDING_2_H
+#define CGAL_FLOAT_SNAP_ROUNDING_2_H
+
+#ifdef DOUBLE_2D_SNAP_VERBOSE
+#include <iostream>
+#endif
+
+#include <CGAL/Surface_sweep_2_algorithms.h>
+
+#include <CGAL/Float_snap_rounding_traits_2.h>
+
+#include <CGAL/intersection_2.h>
+#include <set>
+#include <vector>
+
+#include <CGAL/Named_function_parameters.h>
+
+#include  <CGAL/mutex.h>
+
+namespace CGAL {
+
+namespace internal{
+  /*
+  Scan the vertices from left to right while maintening the y order of the segments.
+  Subdivide the segments if there are too close to a vertex
+  */
+template <class Concurrency_tag=Sequential_tag, class Traits, class PointsRange , class PolylineRange>
+void snap_rounding_scan(PointsRange &pts, PolylineRange &polylines, const Traits &traits){
+  using Point_2   = typename Traits::Point_2;
+  using Segment_2 = typename Traits::Segment_2;
+
+  using Polyline = std::remove_cv_t<typename std::iterator_traits<typename PolylineRange::iterator>::value_type>;
+
+  using Less_xy_2 = typename Traits::Less_xy_2;
+  using Less_y_2 = typename Traits::Less_y_2;
+  using Construct_segment_2 = typename Traits::Construct_segment_2;
+  using Construct_point_at_x_on_segment_2 = typename Traits::Construct_point_at_x_on_segment_2;
+
+  using Compare_squared_distance_2 = typename Traits::Compare_squared_distance_2;
+  using Squared_round_bound_2      = typename Traits::Squared_round_bound_2;
+
+  Less_xy_2 less_xy_2 = traits.less_xy_2_object();
+  Less_y_2 less_y_2 = traits.less_y_2_object();
+  Construct_segment_2 segment_2 = traits.construct_segment_2_object();
+  Construct_point_at_x_on_segment_2 point_at_x = traits.construct_point_at_x_on_segment_2_object();
+
+  Compare_squared_distance_2 csq_dist_2 = traits.compare_squared_distance_2_object();
+  Squared_round_bound_2 round_bound = traits.squared_round_bound_2_object();
+
+  // Precompute round bounds for points
+  std::vector<double> round_bound_pts;
+  round_bound_pts.reserve(pts.size());
+  for(std::size_t i=0; i<pts.size(); ++i)
+    round_bound_pts.emplace_back(round_bound(pts[i]));
+
+  enum class EVENT_TYPE{INSERT, MIDDLE, REMOVE};
+  struct Event{
+    Event(size_t pi_, size_t li_, EVENT_TYPE t_):pi(pi_),li(li_),type(t_){}
+    size_t pi;
+    size_t li;
+    EVENT_TYPE type;
+  };
+
+  auto event_comparator=[&](Event a, Event b){
+    const Point_2 &pa=pts[a.pi];
+    const Point_2 &pb=pts[b.pi];
+
+    // We sort event lexicographically
+    if(a.pi!=b.pi)
+      return less_xy_2(pa,pb);
+
+    // If two events use the same point, we place remove event first for a smaller y_order range
+    if(a.type != b.type)
+      return a.type > b.type;
+    // Only to have a complete order
+    return a.li < b.li;
+  };
+
+  auto pi_below_li=[&](size_t li, std::pair<size_t, size_t> pair_pi_li){
+    std::size_t ind_pi_support = pair_pi_li.second;
+    std::size_t pi = pair_pi_li.first;
+    if(ind_pi_support==li) return false;
+
+    const Polyline &pl=polylines[li];
+
+    Orientation ori=COLLINEAR;
+    if(pl.front()!=pi && pl.back()!=pi)
+      ori=orientation(pts[pl.front()], pts[pl.back()], pts[pi]);
+    // When collinear, pick an other vertex of the line of pi to decide
+    if(ori==COLLINEAR){
+      const Polyline &pi_support = polylines[ind_pi_support];
+      std::size_t other_point = (pi_support.front()!=pi) ? pi_support.front() : pi_support.back();
+      ori=orientation(pts[pl.front()], pts[pl.back()], pts[other_point]);
+      assert(ori!=COLLINEAR);
+    }
+    return ori==RIGHT_TURN;
+  };
+
+#ifdef DOUBLE_2D_SNAP_FULL_VERBOSE
+  std::cout << "Input points" << std::endl;
+  std::size_t i=0;
+  for(const Point_2 &p: pts)
+    std::cout << i++ << ": " << p << std::endl;
+  std::cout << std::endl;
+  std::cout << "Input polylines" << std::endl;
+  i=0;
+  for(const Polyline &pl: polylines){
+    std::cout << i++ << ":";
+    for(std::size_t pi: pl)
+      std::cout << " " << pi;
+    std::cout << std::endl;
+  }
+  std::cout << std::endl;
+#endif
+
+#ifdef DOUBLE_2D_SNAP_VERBOSE
+  std::cout << "Create the Event Queue" << std::endl;
+#endif
+
+  // Input polylines are supposed going from left to right
+  std::vector<Event> event_queue;
+  event_queue.reserve(2*polylines.size());
+  for(size_t li=0; li<polylines.size(); ++li){
+    Polyline &pl=polylines[li];
+    event_queue.emplace_back(pl.front(), li, EVENT_TYPE::INSERT);
+    for(size_t pi=1; pi<pl.size()-1; ++pi){
+      event_queue.emplace_back(pl[pi], li, EVENT_TYPE::MIDDLE);
+      assert(less_xy_2(pts[pl[pi-1]], pts[pl[pi]]));
+    }
+    event_queue.emplace_back(pl.back(), li, EVENT_TYPE::REMOVE);
+    assert(less_xy_2(pts[pl.front()], pts[pl.back()]));
+  }
+#ifdef DOUBLE_2D_SNAP_VERBOSE
+  std::cout << "Sort the event queue along x direction" << std::endl;
+#endif
+  std::sort(event_queue.begin(), event_queue.end(), event_comparator);
+#ifdef DOUBLE_2D_SNAP_VERBOSE
+  std::cout << "Goes through Event" << std::endl;
+#endif
+  // TODO Vector is suboptimal for arbitrary insertion (Actually negligeable in the running time)
+  // Track order of the lines along y along the events
+  std::vector< size_t > y_order;
+  y_order.push_back(event_queue[0].li);
+  for(size_t i=1; i<event_queue.size();++i)
+  {
+    Event event = event_queue[i];
+#ifdef DOUBLE_2D_SNAP_FULL_VERBOSE
+    std::cout << ((event.type==EVENT_TYPE::INSERT)?"Insert": "Remove") << " Event at point " << event.pi << " on line "<< event.li << std::endl;
+#endif
+    // Find the position of the event along the y direction
+    std::vector<size_t>::iterator pos_it = y_order.begin();
+    if(!y_order.empty())
+      pos_it = std::lower_bound(y_order.begin(), y_order.end(), std::pair<size_t, size_t>(event.pi, event.li), pi_below_li);
+    if(event.type==EVENT_TYPE::REMOVE){
+      assert(*pos_it==event.li);
+      pos_it = y_order.erase(pos_it);
+    }
+
+    if(!y_order.empty()){
+      auto close_event=[&](std::size_t pi, std::size_t li){
+        if((pi==polylines[li].front()) || (pi==polylines[li].back()))
+          return true;
+
+        const Point_2 &p = pts[pi];
+        Polyline &pl=polylines[li];
+        Segment_2 seg=segment_2(pts[pl.front()], pts[pl.back()]);
+
+        // (A+B)^2 <= 4*max(A^2,B^2)
+        double bound=round_bound_pts[pi];
+        bound = (std::max)(bound, round_bound_pts[pl.front()]);
+        bound = (std::max)(bound, round_bound_pts[pl.back()]);
+        bound*=4;
+
+        if(possibly(csq_dist_2(p, seg, bound)!=CGAL::LARGER))
+        {
+          if constexpr (std::is_same_v<Exact_predicates_exact_constructions_kernel, typename Traits::Exact_type>)
+          {
+            internal::Evaluate<typename Traits::FT> evaluate;
+            // We refine the pts to reduce the rounding shift and check again
+            evaluate(pts[pi]);
+            evaluate(pts[pl[pl.size()-2]]);
+            // The two following call of exact act on seg variables since they appear in its DAG
+            evaluate(pts[pl.back()]);
+            evaluate(pts[pl.front()]);
+            // Update the bounds
+            round_bound_pts[pi]=round_bound(pts[pi]);
+            round_bound_pts[pl[pl.size()-2]]=round_bound(pts[pl[pl.size()-2]]);
+            round_bound_pts[pl.back()]=round_bound(pts[pl.back()]);
+            round_bound_pts[pl.front()]=round_bound(pts[pl.front()]);
+            bound=round_bound_pts[pi];
+            bound = (std::max)(bound, round_bound_pts[pl[pl.size()-2]]);
+            bound = (std::max)(bound, round_bound_pts[pl.back()]);
+            bound*=4;
+
+            // Check if the point and the segment are still too closed for a safe rounding
+            if(csq_dist_2(p, seg, bound)==CGAL::LARGER)
+              return false;
+          }
+
+          // Check if segment was not already subdivided by another point
+          for(std::size_t i: pl)
+            if(pts[i].x()==pts[pi].x())
+              return false;
+
+          // Create a point on seg at the same x coordinate than p
+          pts.push_back(point_at_x(seg, pts[event.pi].x()));
+          std::size_t new_pi=pts.size()-1;
+          round_bound_pts.emplace_back(round_bound(pts[new_pi]));
+          // We insert it on pl before the last vertex
+          pl.insert(pl.end()-1, new_pi);
+#ifdef DOUBLE_2D_SNAP_FULL_VERBOSE
+          std::cout << "Create point " << new_pi << " on " << li << " due to proximity with " << pi << "_____________________________" << std::endl;
+          std::cout << new_pi <<": " << pts[new_pi] << std::endl;
+          std::cout << li << ":";
+          for(std::size_t i: pl)
+            std::cout << " " << i;
+          std::cout << std::endl;
+#endif
+          return true;
+        }
+        return false;
+      };
+
+      if(event.pi != event_queue[i-1].pi)
+      {
+        // Look above segments and creates a point if there too close for a safe rounding
+        auto above = pos_it;
+        size_t pi=event.pi;
+        while(above!=y_order.end() && close_event(pi,*above)){
+          if((pi!=polylines[*above].front()) && (pi!=polylines[*above].back()))
+            pi=pts.size()-1;
+          ++above;
+        }
+
+        // same with below segments
+        auto below = pos_it;
+        pi=event.pi;
+        if(below!=y_order.begin()){
+          --below;
+          while(close_event(pi,*(below))){
+            if(below==y_order.begin())
+              break;
+            if((pi!=polylines[*below].front()) && (pi!=polylines[*below].back()))
+              pi=pts.size()-1;
+            --below;
+          }
+        }
+      }
+    }
+
+    if(event.type==EVENT_TYPE::INSERT){
+      y_order.insert(pos_it, event.li);
+    }
+  }
+
+  // We sort the point by y to ensure the order of the point is preserved
+  // Round value of y coordinates are indirectly modified in the case of a filter failure
+  auto sort_pi=[&](std::size_t i, std::size_t j){
+    return less_y_2(pts[i],pts[j]);
+  };
+  std::vector< std::size_t > indices(pts.size(),0);
+  std::iota(indices.begin(),indices.end(),0);
+  std::sort(indices.begin(),indices.end(),sort_pi);
+
+#ifdef DOUBLE_2D_SNAP_FULL_VERBOSE
+  std::cout << "Output points" << std::endl;
+  i=0;
+  for(const Point_2 &p: pts)
+    std::cout << i++ << ": " << p << std::endl;
+  std::cout << std::endl;
+  std::cout << "Output polylines" << std::endl;
+  i=0;
+  for(const Polyline &pl: polylines){
+    std::cout << i++ << ":";
+    for(std::size_t pi: pl)
+      std::cout << " " << pi;
+    std::cout << std::endl;
+  }
+  std::cout << std::endl;
+#endif
+}
+
+template <class Concurrency_tag=Sequential_tag, class Traits, class PointsRange , class PolylineRange>
+void merge_duplicate_points_in_polylines(PointsRange &pts, PolylineRange &polylines, const Traits &traits)
+{
+  using Point_2 = typename Traits::Point_2;
+  using Polyline = std::remove_cv_t<typename std::iterator_traits<typename PolylineRange::iterator>::value_type>;
+
+  using Less_xy_2 = typename Traits::Less_xy_2;
+  using Equal_2 = typename Traits::Equal_2;
+  Equal_2 equal = traits.equal_2_object();
+
+  auto Less_indexes_xy_2=[&](std::size_t i, std::size_t j){
+    return Less_xy_2()(pts[i], pts[j]);
+  };
+
+  std::vector< std::size_t > unique_points(pts.size());
+  std::iota(unique_points.begin(), unique_points.end(), 0);
+  std::sort(unique_points.begin(), unique_points.end(), Less_indexes_xy_2);
+  std::vector<Point_2> new_pts;
+  std::vector<std::size_t> old_to_new_index(pts.size());
+  for(std::size_t i=0; i!=pts.size(); ++i){
+    if(i==0 || !equal(pts[unique_points[i]],pts[unique_points[i-1]]))
+      new_pts.push_back(pts[unique_points[i]]);
+    old_to_new_index[unique_points[i]]=new_pts.size()-1;
+  }
+
+  std::swap(pts, new_pts);
+  for (Polyline& polyline : polylines) {
+    std::vector<std::size_t> updated_polyline;
+    for (std::size_t i=0; i<polyline.size(); ++i) {
+        std::size_t new_pi=old_to_new_index[polyline[i]];
+        if(i==0 || (new_pi!=updated_polyline[updated_polyline.size()-1]))
+          updated_polyline.push_back(new_pi);
+        assert(new_pi<pts.size());
+        assert(pts[new_pi]==new_pts[polyline[i]]);
+    }
+    std::swap(polyline, updated_polyline);
+  }
+}
+
+// Some points may have collapsed on a vertical segment, we subdivide these vertical segments accordingly
+template <class Concurrency_tag=Sequential_tag, class Traits, class PointsRange , class PolylineRange>
+void snap_post_process(PointsRange &pts, PolylineRange &polylines, const Traits &traits)
+{
+  using Polyline = std::remove_cv_t<typename std::iterator_traits<typename PolylineRange::iterator>::value_type>;
+  using Less_xy_2 = typename Traits::Less_xy_2;
+
+  Less_xy_2 less_xy_2 = traits.less_xy_2_object();
+
+  auto Less_indexes_xy_2=[&](std::size_t i, std::size_t j){
+    return less_xy_2(pts[i], pts[j]);
+  };
+
+  std::set<std::size_t, decltype(Less_indexes_xy_2)> p_sort_by_x(Less_indexes_xy_2);
+  for(std::size_t i=0; i!=pts.size(); ++i)
+    p_sort_by_x.insert(i);
+
+  using Iterator_set_x = typename std::set<std::size_t, decltype(Less_indexes_xy_2)>::iterator;
+
+  for(Polyline &poly: polylines){
+    std::vector<std::size_t> updated_polyline;
+    updated_polyline.push_back(poly.front());
+    for(std::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(Less_indexes_xy_2(poly[i-1],poly[i])){
+          start=p_sort_by_x.upper_bound(poly[i-1]);
+          end=p_sort_by_x.lower_bound(poly[i]);
+          for(auto it=start; it!=end; ++it){
+            updated_polyline.push_back(*it);
+          }
+        } else {
+          start=p_sort_by_x.upper_bound(poly[i]);
+          end=p_sort_by_x.lower_bound(poly[i-1]);
+          for(auto it=end; it!=start;){
+            --it;
+            updated_polyline.push_back(*it);
+          }
+        }
+      }
+      updated_polyline.push_back(poly[i]);
+    }
+    std::swap(poly, updated_polyline);
+  }
+}
+
+template <class Concurrency_tag=Sequential_tag, class Traits, class PointsRange , class PolylineRange>
+void double_snap_rounding_2_disjoint(PointsRange &pts, PolylineRange &polylines, const Traits &traits)
+{
+  using Point_2   = typename Traits::Point_2;
+  using Construct_round_point_2    = typename Traits::Construct_round_point_2;
+
+  Construct_round_point_2 round = traits.construct_round_point_2_object();
+
+  snap_rounding_scan(pts, polylines, traits);
+#ifdef DOUBLE_2D_SNAP_VERBOSE
+  std::cout << "Round" << std::endl;
+#endif
+  for(Point_2 &p: pts)
+    p=round(p);
+#ifdef DOUBLE_2D_SNAP_VERBOSE
+  std::cout << "Merging points that are collapsed together" << std::endl;
+#endif
+  merge_duplicate_points_in_polylines(pts, polylines, traits);
+#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 vertical segments with vertices on them" << std::endl;
+#endif
+  snap_post_process(pts, polylines, traits);
+}
+
+} // end of namespace internal
+
+/**
+* ingroup PkgSnapRounding2Ref
+*
+* Given a range of segments, compute rounded subsegments that are pairwise disjoint in their interiors, based on the input curves.
+* The output is a range of polyline with each polyline corresponding to an input segment.
+*
+* TODO Currently compute_subcurves have no visitor to obtain a polyline per input segment, thus
+* currently the output contain one polyline per ubsegment of the arrangement
+*
+* @tparam Concurrency_tag allows choosing whether the algorithm runs in parallel or sequentially.
+* If `CGAL::Parallel_tag` is specified and CGAL is linked with the Intel TBB library, the algorithm will run in parallel.
+* Otherwise, if `CGAL::Sequential_tag` is specified (the default), the algorithm will run sequentially.
+* @tparam InputIterator iterator over a range of `Segment_2`
+* @tparam OutputContainer inserter over a range of `Polyline`. `Polyline` must be a type that provides a `push_back(Point_2)` function.
+* @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
+* \param begin,end the input segment range
+* \param out the output container
+* \param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+* \cgalNamedParamsBegin
+*   \cgalParamNBegin{concurrency_tag}
+*     \cgalParamDescription{That template parameter enables to choose whether the algorithm is to be run in parallel, if CGAL::Parallel_tag
+*                           is specified and CGAL has been linked with the Intel TBB library, or sequentially, otherwise.}
+*     \cgalParamType{CGAL::Concurrency_tag}
+*     \cgalParamDefault{CGAL::Sequential_tag}
+*   \cgalParamNEnd
+*
+*   \cgalParamNBegin{geom_traits}
+*     \cgalParamDescription{an instance of a geometric traits class}
+*     \cgalParamType{The traits class must respect the concept of `FloatSnapRoundingTraits_2`}
+*     \cgalParamDefault{an instance of `Float_snap_rounding_traits_2`}
+*   \cgalParamNEnd
+*/
+template <class InputIterator , class OutputContainer, class NamedParameters = parameters::Default_named_parameters>
+typename OutputContainer::iterator double_snap_rounding_2(InputIterator  	begin,
+		                                                      InputIterator  	end,
+		                                                      OutputContainer &out,
+                                                          const NamedParameters &np = parameters::default_values())
+{
+  using Concurrency_tag = typename internal_np::Lookup_named_param_def<internal_np::concurrency_tag_t,
+                                                              NamedParameters,
+                                                              Sequential_tag>::type;
+
+  using InputKernel = typename Kernel_traits<std::remove_cv_t<typename std::iterator_traits<InputIterator>::value_type>>::Kernel;
+  using DefaultTraits = Float_snap_rounding_traits_2<InputKernel>;
+  using Traits = typename internal_np::Lookup_named_param_def<internal_np::geom_traits_t,
+                                                              NamedParameters,
+                                                              DefaultTraits>::type;
+
+  using Point_2 = typename Traits::Point_2;
+  using Segment_2 = typename Traits::Segment_2;
+  using I2E = typename Traits::Converter_to_exact;
+  using E2O = typename Traits::Converter_from_exact;
+  using VectorIterator = typename std::vector<Segment_2>::iterator;
+
+  using Polyline = std::remove_cv_t<typename std::iterator_traits<typename OutputContainer::iterator>::value_type>;
+
+  using Less_xy_2 = typename Traits::Less_xy_2;
+
+  using parameters::choose_parameter;
+  using parameters::get_parameter;
+
+  const Traits &traits = choose_parameter(get_parameter(np, internal_np::geom_traits), Traits());
+
+  I2E to_exact=traits.converter_to_exact_object();
+  E2O from_exact=traits.converter_from_exact_object();
+
+  std::vector<Segment_2> convert_input;
+  for(InputIterator it=begin; it!=end; ++it)
+    convert_input.push_back(to_exact(*it));
+  std::vector<Segment_2> segs;
+#ifdef DOUBLE_2D_SNAP_VERBOSE
+  std::cout << "Solved intersections" << std::endl;
+  std::cout << "do intersect? " << do_curves_intersect(convert_input.begin(), convert_input.end()) << std::endl;
+#endif
+  compute_subcurves(convert_input.begin(), convert_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< std::size_t> > polylines;
+
+  // Transform range of the segments in the range of points and polyline of indexes
+  for(VectorIterator it=segs.begin(); it!=segs.end(); ++it)
+  {
+    const Point_2& p1 = it->source();
+    const Point_2& p2 = it->target();
+
+    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(VectorIterator it=segs.begin(); it!=segs.end(); ++it)
+  {
+    std::size_t index1 = point_to_index[it->source()];
+    std::size_t index2 = point_to_index[it->target()];
+    if(Less_xy_2()(it->source(), it->target()))
+      polylines.push_back({index1, index2});
+    else
+      polylines.push_back({index2, index1});
+  }
+
+  // Main algorithm
+  internal::double_snap_rounding_2_disjoint<Concurrency_tag, Traits>(pts, polylines, traits);
+
+#ifdef DOUBLE_2D_SNAP_VERBOSE
+  std::cout << "Build output" << std::endl;
+#endif
+
+  // Output polylines
+  for(auto &poly: polylines){
+    Polyline new_line;
+    for(std::size_t pi: poly)
+      new_line.push_back(from_exact(pts[pi]));
+    out.push_back(new_line);
+  }
+
+  return out.end();
+}
+
+/**
+* ingroup PkgSnapRounding2Ref
+*
+* Given a range of segments, compute rounded subsegments that are pairwise disjoint in their interior, as induced by the input curves.
+*
+* @tparam Concurrency_tag That template parameter enables to choose whether the algorithm is to be run in
+* parallel, if CGAL::Parallel_tag is specified and CGAL has been linked with the Intel TBB library, or sequentially, if CGAL::Sequential_tag - the default value - is specified.
+* @tparam InputIterator iterator of a segment range
+* @tparam OutputContainer inserter of a segment range
+* @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
+* \param begin,end the input segment range
+* \param out the output inserter
+* \param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+* \cgalNamedParamsBegin
+*   \cgalParamNBegin{concurrency_tag}
+*     \cgalParamDescription{That template parameter enables to choose whether the algorithm is to be run in parallel, if CGAL::Parallel_tag
+*                           is specified and CGAL has been linked with the Intel TBB library, or sequentially, otherwise.}
+*     \cgalParamType{CGAL::Concurrency_tag}
+*     \cgalParamDefault{CGAL::Sequential_tag}
+*   \cgalParamNEnd
+*
+*   \cgalParamNBegin{geom_traits}
+*     \cgalParamDescription{an instance of a geometric traits class}
+*     \cgalParamType{The traits class must respect the concept of `FloatSnapRoundingTraits_2`}
+*     \cgalParamDefault{an instance of `Float_snap_rounding_traits_2`}
+*   \cgalParamNEnd
+*/
+template <class InputIterator , class OutputIterator, class NamedParameters = parameters::Default_named_parameters>
+OutputIterator compute_snapped_subcurves_2(InputIterator  	 begin,
+		                                       InputIterator  	 end,
+		                                       OutputIterator    out,
+                                           const NamedParameters &np = parameters::default_values())
+{
+  using Concurrency_tag = typename internal_np::Lookup_named_param_def<internal_np::concurrency_tag_t,
+                                                              NamedParameters,
+                                                              Sequential_tag>::type;
+
+  using InputKernel = typename Kernel_traits<std::remove_cv_t<typename std::iterator_traits<InputIterator>::value_type>>::Kernel;
+  using DefaultTraits = Float_snap_rounding_traits_2<InputKernel>;
+  using Traits = typename internal_np::Lookup_named_param_def<internal_np::geom_traits_t,
+                                                              NamedParameters,
+                                                              DefaultTraits>::type;
+
+  using Point_2 = typename Traits::Point_2;
+  using Segment_2 = typename Traits::Segment_2;
+  using I2E = typename Traits::Converter_to_exact;
+  using E2O = typename Traits::Converter_from_exact;
+
+  using SegmentRange = std::vector<Segment_2>;
+  using SegmentRangeIterator = typename SegmentRange::iterator;
+
+  using Less_xy_2 = typename Traits::Less_xy_2;
+  using Construct_segment_2 = typename Traits::Construct_segment_2;
+
+  using parameters::choose_parameter;
+  using parameters::get_parameter;
+
+  const Traits &traits = choose_parameter(get_parameter(np, internal_np::geom_traits), Traits());
+
+  I2E to_exact=traits.converter_to_exact_object();
+  E2O from_exact=traits.converter_from_exact_object();
+
+  Less_xy_2 less_xy_2 = traits.less_xy_2_object();
+  Construct_segment_2 segment_2 = traits.construct_segment_2_object();
+
+  std::vector<Segment_2> convert_input;
+  for(InputIterator it=begin; it!=end; ++it)
+    if(it->source()!=it->target())
+      convert_input.push_back(to_exact(*it));
+  std::vector<Segment_2> segs;
+#ifdef DOUBLE_2D_SNAP_VERBOSE
+  std::cout << "Solved intersections" << std::endl;
+  std::cout << "do intersect? " << do_curves_intersect(convert_input.begin(), convert_input.end()) << std::endl;
+#endif
+  compute_subcurves(convert_input.begin(), convert_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, Less_xy_2> unique_point_set(less_xy_2);
+  std::map<Point_2, std::size_t> point_to_index;
+  std::vector<Point_2> pts;
+  std::vector< std::vector< std::size_t> > polylines;
+
+  // Transform range of the segments in the range of points and polyline of indexes
+  for(SegmentRangeIterator it=segs.begin(); it!=segs.end(); ++it)
+  {
+    const Point_2& p1 = it->source();
+    const Point_2& p2 = it->target();
+
+    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(SegmentRangeIterator it=segs.begin(); it!=segs.end(); ++it)
+  {
+    std::size_t index1 = point_to_index[it->source()];
+    std::size_t index2 = point_to_index[it->target()];
+    if(less_xy_2(it->source(), it->target()))
+      polylines.push_back({index1, index2});
+    else
+      polylines.push_back({index2, index1});
+  }
+
+
+  // Main algorithm
+  internal::double_snap_rounding_2_disjoint<Concurrency_tag>(pts, polylines, traits);
+
+#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<std::size_t,std::size_t> > set_out_segs;
+  for(auto &poly: polylines)
+    for(std::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){
+    *out++=from_exact(segment_2(pts[pair.first], pts[pair.second]));
+    assert(pts[pair.first]!=pts[pair.second]);
+  }
+
+  return out;
+}
+
+/**
+* ingroup PkgSnapRounding2Ref
+*
+* Given a range of `Polygon_2`, compute rounded polygons such that their segments are either equal either disjoint in their interior, as induced by the input polygons.
+* The polygons are intended to be non-intersecting, unless the named parameter `compute_intersection` is set to `true`.
+* Any polygon is guarantee to remain a Polygon in the output but may present pinched section or/and common vertices or segments with
+* other polygons.
+*
+* @tparam InputIterator iterator of a CGAL::Polygon_2 range
+* @tparam OutputContainer inserter of a CGAL::Polygon_2 range
+* @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
+* \param begin,end the input polygon range
+* \param out the output inserter
+* \param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+* \cgalNamedParamsBegin
+*   \cgalParamNBegin{concurrency_tag}
+*     \cgalParamDescription{That template parameter enables to choose whether the algorithm is to be run in parallel, if CGAL::Parallel_tag
+*                           is specified and CGAL has been linked with the Intel TBB library, or sequentially, otherwise.}
+*     \cgalParamType{CGAL::Concurrency_tag}
+*     \cgalParamDefault{CGAL::Sequential_tag}
+*   \cgalParamNEnd
+*
+*   \cgalParamNBegin{compute_intersection}
+*     \cgalParamDescription{Enable intersection computation between the polygons before performing snapping.}
+*     \cgalParamType{boolean}
+*     \cgalParamDefault{false}
+*   \cgalParamNEnd
+*
+*   \cgalParamNBegin{geom_traits}
+*     \cgalParamDescription{an instance of a geometric traits class}
+*     \cgalParamType{The traits class must respect the concept of `FloatSnapRoundingTraits_2`}
+*     \cgalParamDefault{an instance of `Float_snap_rounding_traits_2`}
+*   \cgalParamNEnd
+* @warning The convex property of the polygons is not necessarly preserved
+*/
+template <class InputIterator, class OutputIterator, class NamedParameters = parameters::Default_named_parameters>
+void compute_snapped_polygons_2(InputIterator begin,
+                                InputIterator end,
+                                OutputIterator out,
+                                const NamedParameters &np = parameters::default_values())
+{
+  using Concurrency_tag = typename internal_np::Lookup_named_param_def<internal_np::concurrency_tag_t,
+                                                              NamedParameters,
+                                                              Sequential_tag>::type;
+
+  using Polygon_2 = typename std::iterator_traits<InputIterator>::value_type;
+  using InputKernel = typename Kernel_traits<typename Polygon_2::Point_2>::Kernel;
+  using DefaultTraits = Float_snap_rounding_traits_2<InputKernel>;
+  using Traits = typename internal_np::Lookup_named_param_def<internal_np::geom_traits_t,
+                                                              NamedParameters,
+                                                              DefaultTraits>::type;
+
+  using Less_xy_2 = typename Traits::Less_xy_2;
+
+  using Point_2 = typename Traits::Point_2;
+  using I2E = typename Traits::Converter_to_exact;
+  using E2O = typename Traits::Converter_from_exact;
+
+  using parameters::choose_parameter;
+  using parameters::get_parameter;
+
+  Polygon_2 P=*begin;
+
+  const Traits &traits = choose_parameter(get_parameter(np, internal_np::geom_traits), Traits());
+  const bool compute_intersections = choose_parameter(get_parameter(np, internal_np::compute_intersection), false);
+
+  Less_xy_2 less_xy_2 = traits.less_xy_2_object();
+
+  I2E to_exact=traits.converter_to_exact_object();
+  E2O from_exact=traits.converter_from_exact_object();
+
+#ifdef DOUBLE_2D_SNAP_VERBOSE
+  std::cout << "Change format to range of points and indexes" << std::endl;
+#endif
+  std::vector<Point_2> pts;
+  std::vector< std::vector< std::size_t> > polylines;
+  std::vector<size_t> polygon_index;
+
+  if(compute_intersections){
+    // TODO Need to compute_subcurves that output polylines to track the polygons
+    assert(0);
+  } else {
+    // Index of the segments that introduced a new polygon
+    polygon_index.reserve(std::distance(begin, end));
+    for(InputIterator it=begin; it!=end; ++it){
+      size_t index_start=polylines.size();
+      polygon_index.push_back(polylines.size());
+      for(const typename Polygon_2::Point_2 &p: it->vertices())
+        pts.push_back(to_exact(p));
+      for(size_t i=0; i<P.size()-1; ++i)
+        if(less_xy_2(pts[i], pts[i+1]))
+          polylines.push_back({i, i+1});
+        else
+          polylines.push_back({i+1, i});
+      polylines.push_back({pts.size()-1,index_start});
+      assert(pts.size()==polylines.size());
+    }
+    polygon_index.push_back(polylines.size());
+  }
+
+  // Main algorithm
+  internal::double_snap_rounding_2_disjoint<Concurrency_tag>(pts, polylines, traits);
+
+#ifdef DOUBLE_2D_SNAP_VERBOSE
+  std::cout << "Build output" << std::endl;
+#endif
+
+  // Output a range of segments while removing duplicate ones
+  for(std::size_t input_ind=0; input_ind<std::size_t(std::distance(begin,end)); ++input_ind){
+    for(std::size_t pl_ind=polygon_index[input_ind]; pl_ind<polygon_index[input_ind+1]; ++pl_ind){
+      std::vector<std::size_t> &poly = polylines[pl_ind];
+      Polygon_2 P;
+      for(std::size_t i=1; i<poly.size(); ++i)
+        P.push_back(from_exact(pts[poly[i]]));
+    }
+    *out++=P;
+  }
+}
+
+
+/**
+* ingroup
+*
+* Given a Polygon_2, compute rounded segments that are pairwise disjoint in their interior, as induced by the input polygon.
+* The output is guarantee to be a Polygon but may present pinched section.
+*
+* \tparam Polygon_2 model of CGAL::Polygon_2
+* \tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
+* \param P the input polygon
+* \param out the output polygon
+* \param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+* \cgalNamedParamsBegin
+*   \cgalParamNBegin{concurrency_tag}
+*     \cgalParamDescription{That template parameter enables to choose whether the algorithm is to be run in parallel, if CGAL::Parallel_tag
+*                           is specified and CGAL has been linked with the Intel TBB library, or sequentially, otherwise.}
+*     \cgalParamType{CGAL::Concurrency_tag}
+*     \cgalParamDefault{CGAL::Sequential_tag}
+*   \cgalParamNEnd
+*
+*   \cgalParamNBegin{geom_traits}
+*     \cgalParamDescription{a multiplier of the error value for boundary edges to preserve the boundaries}
+*     \cgalParamType{double}
+*     \cgalParamDefault{100}
+*   \cgalParamNEnd
+* @warning The convex property is not necessarly preserved
+*/
+template <class Polygon_2, class NamedParameters = parameters::Default_named_parameters>
+void compute_snapped_polygon_2(const Polygon_2 &P, Polygon_2 &out, const NamedParameters &np = parameters::default_values())
+{
+  std::array<Polygon_2, 1> vec({P});
+  std::vector<Polygon_2> out_vec;
+  compute_snapped_polygons_2(vec.begin(), vec.end(), std::back_inserter(out_vec), np);
+  out = out_vec[0];
+}
+
+
+
+} //namespace CGAL
+
+#endif
diff --git a/Snap_rounding_2/include/CGAL/Float_snap_rounding_traits_2.h b/Snap_rounding_2/include/CGAL/Float_snap_rounding_traits_2.h
new file mode 100644
index 00000000000..1d61398de05
--- /dev/null
+++ b/Snap_rounding_2/include/CGAL/Float_snap_rounding_traits_2.h
@@ -0,0 +1,116 @@
+// Copyright (c) 2025 Geometry Factory.
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+//
+// author(s)     : Léo Valque
+
+#ifndef CGAL_FLOAT_SNAP_ROUNDING_TRAITS_2_H
+#define CGAL_FLOAT_SNAP_ROUNDING_TRAITS_2_H
+
+#include <CGAL/Arr_segment_traits_2.h>
+
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Cartesian_converter.h>
+
+#include <CGAL/Named_function_parameters.h>
+
+#include  <CGAL/mutex.h>
+
+namespace CGAL {
+
+/*!
+\ingroup PkgSnapRounding2Ref
+
+The class `Float_snap_rounding_traits_2<Kernel>` is a model of the
+`FloatSnapRoundingTraits_2` concept, and is the only traits class supplied
+with the package.
+This class should be instantiated with a kernel with points and segments that are convertible to
+`Exact_predicates_exact_constructions_kernel`, the user must specified to used another exact geometric kernel that conforms
+to the \cgal kernel-concept, such as `Cartesian<Gmpq>`.
+
+\cgalModels{SnapRoundingTraits_2}
+
+*/
+template<typename Input_Kernel, typename Exact_Kernel = Exact_predicates_exact_constructions_kernel>
+struct Float_snap_rounding_traits_2: Arr_segment_traits_2<Exact_Kernel>{
+  typedef Arr_segment_traits_2<Exact_Kernel> Base;
+
+  typedef Exact_Kernel Exact_type;
+
+  typedef typename Base::FT        FT;
+  typedef typename Base::Point_2   Point_2;
+  typedef typename Base::Segment_2 Segment_2;
+
+  typedef typename Base::Less_xy_2 Less_xy_2;
+  typedef typename Base::Less_y_2  Less_y_2;
+  typedef typename Base::Equal_2   Equal_2;
+
+  typedef typename Base::Construct_point_2   Construct_point_2;
+  typedef typename Base::Construct_source_2  Construct_source_2;
+  typedef typename Base::Construct_target_2  Construct_target_2;
+  typedef typename Base::Construct_segment_2 Construct_segment_2;
+
+  typedef Cartesian_converter<Input_Kernel, Exact_Kernel> Converter_to_exact;
+  typedef Cartesian_converter<Exact_Kernel, Input_Kernel> Converter_from_exact;
+
+  // Return an upperbound of the squared distance between a point and its rounded value
+  struct Squared_round_bound_2{
+    double operator()(const FT &x) const{
+      double b=std::nextafter(to_interval(x).second - to_interval(x).first, std::numeric_limits<double>::infinity());
+      return b*b;
+    }
+    double operator()(const Point_2 &p) const{
+      return (*this)(p.x())+(*this)(p.y());
+    }
+    double operator()(const Segment_2 &seg) const{
+      return (std::max)( (*this)(Base().construct_source_2_object()(seg)),
+                         (*this)(Base().construct_target_2_object()(seg)));
+    }
+  };
+
+  struct Construct_round_point_2{
+    double operator()(const FT &x) const{
+      return to_double(x);
+    }
+    Point_2 operator()(const Point_2 &p) const{
+      return Point_2((*this)(p.x()),(*this)(p.y()));
+    }
+  };
+
+  struct Construct_point_at_x_on_segment_2{
+    Point_2 operator()(const Segment_2 &seg, const FT &x) const{
+      FT y= (seg.supporting_line().y_at_x(x));
+      return Base().construct_point_2_object()(x, y);
+    }
+  };
+
+  Converter_to_exact converter_to_exact_object() const{
+    return Converter_to_exact();
+  }
+
+  Converter_from_exact converter_from_exact_object() const{
+    return Converter_from_exact();
+  }
+
+  Squared_round_bound_2 squared_round_bound_2_object() const{
+    return Squared_round_bound_2();
+  }
+
+  Construct_round_point_2 construct_round_point_2_object() const{
+    return Construct_round_point_2();
+  }
+
+  Construct_point_at_x_on_segment_2 construct_point_at_x_on_segment_2_object() const{
+    return Construct_point_at_x_on_segment_2();
+  }
+};
+
+} //namespace CGAL
+
+#endif
diff --git a/Snap_rounding_2/include/CGAL/Snap_rounding_traits_2.h b/Snap_rounding_2/include/CGAL/Snap_rounding_traits_2.h
index 517c56907a3..c7d6d16b7cb 100644
--- a/Snap_rounding_2/include/CGAL/Snap_rounding_traits_2.h
+++ b/Snap_rounding_2/include/CGAL/Snap_rounding_traits_2.h
@@ -59,6 +59,10 @@ public:
   public:
     void operator()(const Point_2& p, NT pixel_size, NT &x, NT &y)
     {
+      // TODO a proper and efficient double_ceil was done for snap_rounding_polygon_soup
+#if 0
+      x = double_ceil((x / pixel_size) - 0.5) * pixel_size;
+#else
       NT x_tmp = p.x() / pixel_size;
       NT y_tmp = p.y() / pixel_size;
 
@@ -81,6 +85,7 @@ public:
 
       x = x_floor * pixel_size + pixel_size / NT(2.0);
       y = y_floor * pixel_size + pixel_size / NT(2.0);
+#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..94bc16ced5c 100644
--- a/Snap_rounding_2/test/Snap_rounding_2/CMakeLists.txt
+++ b/Snap_rounding_2/test/Snap_rounding_2/CMakeLists.txt
@@ -7,6 +7,8 @@ project(Snap_rounding_2_Tests)
 find_package(CGAL REQUIRED)
 
 create_single_source_cgal_program(test_snap_rounding_2.cpp NO_TESTING)
+create_single_source_cgal_program(test_float_snap_rounding_2.cpp)
+create_single_source_cgal_program(test_float_snap_rounding_2_concept.cpp)
 
 function(add_Snap_rounding_tests name)
   set(data_dir "data")
diff --git a/Snap_rounding_2/test/Snap_rounding_2/test_float_snap_rounding_2.cpp b/Snap_rounding_2/test/Snap_rounding_2/test_float_snap_rounding_2.cpp
new file mode 100644
index 00000000000..421204fbde0
--- /dev/null
+++ b/Snap_rounding_2/test/Snap_rounding_2/test_float_snap_rounding_2.cpp
@@ -0,0 +1,331 @@
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/Float_snap_rounding_traits_2.h>
+#include <CGAL/Float_snap_rounding_2.h>
+#include <CGAL/Arr_segment_traits_2.h>
+#include <CGAL/Surface_sweep_2_algorithms.h>
+
+#include <CGAL/Boolean_set_operations_2.h>
+
+#include <CGAL/Random.h>
+
+typedef CGAL::Exact_predicates_exact_constructions_kernel       Kernel;
+typedef CGAL::Float_snap_rounding_traits_2<Kernel>              Traits_2;
+typedef Kernel::Segment_2                                       Segment_2;
+typedef Kernel::Point_2                                         Point_2;
+typedef Kernel::Vector_2                                        Vector_2;
+typedef std::vector<Point_2 >                                   Polyline_2;
+typedef std::vector<Polyline_2>                                 Polyline_range_2;
+typedef Kernel::FT                                              FT;
+
+typedef CGAL::Arr_segment_traits_2<Kernel>                      Arr_segment_traits_2;
+typedef Arr_segment_traits_2::Curve_2                           Curve_2;
+
+typedef CGAL::Polygon_2<Kernel>                                 Polygon_2;
+typedef CGAL::Polygon_with_holes_2<Kernel>                      Polygon_with_holes_2;
+typedef std::vector<Polygon_with_holes_2>                       Pwh_vec_2;
+
+#ifdef BENCH_AND_VERBOSE_FLOAT_SNAP_ROUNDING_2
+#include <CGAL/Real_timer.h>
+#include <CGAL/Cartesian_converter.h>
+#include <CGAL/Simple_cartesian.h>
+typedef CGAL::Simple_cartesian<double>                          Naive;
+typedef CGAL::Cartesian_converter<Kernel,Naive>                 EK_to_IK;
+typedef CGAL::Cartesian_converter<Naive, Kernel>                IK_to_EK;
+#endif
+
+#ifdef COMPARE_WITH_INTEGER_SNAP_ROUNDING_2
+#include <CGAL/Snap_rounding_traits_2.h>
+#include <CGAL/Snap_rounding_2.h>
+typedef CGAL::Snap_rounding_traits_2<Kernel>                     SnapTraits;
+#endif
+
+//Biggest double with ulp smaller than an integer
+const double maxFloat = std::pow(2,23);
+const double maxDouble = std::pow(2,52);
+
+Segment_2 scale_segment(const Segment_2 &s){
+  auto scale_value=[](FT x){
+    return maxDouble * CGAL::to_double(x);
+  };
+  return Segment_2(Point_2(scale_value(s.start().x()), scale_value(s.start().y())),
+                   Point_2(scale_value(s.end().x()), scale_value(s.end().y())));
+}
+
+Point_2 random_point(CGAL::Random& r, double m = -1, double M = 1)
+{
+  return Point_2(FT(r.get_double(m, M)), FT(r.get_double(m, M)));
+}
+
+Point_2 random_noise_point(CGAL::Random& r, double m, double M, Vector_2 v)
+{
+  return Point_2(FT(r.get_double(m, M) + CGAL::to_double(v.x())), FT(r.get_double(m, M)  + CGAL::to_double(v.y())));
+}
+
+#ifdef BENCH_AND_VERBOSE_FLOAT_SNAP_ROUNDING_2
+template<class Iterator, class OutRange>
+void compute_subcurves_and_naive_snap(Iterator begin, Iterator end, OutRange &out){
+  EK_to_IK to_inexact;
+  IK_to_EK to_exact;
+  std::vector<Curve_2> segs;
+  for(Iterator it=begin; it!=end; ++it)
+    if(it->source()!=it->target())
+      segs.emplace_back(*it);
+  CGAL::compute_subcurves(segs.begin(), segs.end(), std::back_inserter(out));
+  //Naive snap
+  OutRange out2;
+  for(Segment_2 &seg: out){
+    Segment_2 s = to_exact(to_inexact(seg));
+    if(!s.is_degenerate())
+      out2.emplace_back(s);
+  }
+  std::swap(out,out2);
+}
+#endif
+
+void test(const std::vector<Segment_2> &segs){
+  std::vector<Segment_2> out;
+#ifdef BENCH_AND_VERBOSE_FLOAT_SNAP_ROUNDING_2
+  CGAL::Real_timer t;
+  std::cout << "Input size: " << segs.size() << std::endl;
+  t.start();
+  compute_subcurves_and_naive_snap(segs.begin(), segs.end(), out);
+  t.stop();
+  std::cout << "Naive snap size: " << out.size() << " ,running time: " << t.time() << " and output do intersect: " << CGAL::do_curves_intersect(out.begin(), out.end()) <<std::endl;
+  out.clear();
+  t.reset();
+  t.start();
+#endif
+  CGAL::compute_snapped_subcurves_2(segs.begin(), segs.end(), std::back_inserter(out));
+  std::vector<Segment_2> out2(out.begin(), out.end());
+#ifdef BENCH_AND_VERBOSE_FLOAT_SNAP_ROUNDING_2
+  t.stop();
+  std::cout << "Formal snap size: " << out.size() << " ,running time: " << t.time() << std::endl;
+#endif
+  assert(!CGAL::do_curves_intersect(out.begin(), out.end()));
+#ifdef COMPARE_WITH_INTEGER_SNAP_ROUNDING_2
+  Polyline_range_2 output_list;
+  t.reset();
+  t.start();
+  CGAL::snap_rounding_2<SnapTraits>(segs.begin(), segs.end(), std::back_inserter(output_list), 1./maxDouble);
+  t.stop();
+  std::cout << "Running time with integer 2D snap (scaled at 10^15): " << t.time() << std::endl;
+#endif
+#ifdef BENCH_AND_VERBOSE_FLOAT_SNAP_ROUNDING_2
+  std::cout << "\n";
+#endif
+}
+
+void test_fully_random(CGAL::Random &r, size_t nb_segments){
+#ifdef BENCH_AND_VERBOSE_FLOAT_SNAP_ROUNDING_2
+  std::cout << "Test fully random" << std::endl;
+#endif
+  std::vector<Segment_2> segs;
+  for(size_t i=0; i!=nb_segments; ++i)
+    segs.emplace_back(random_point(r), random_point(r));
+
+  test(segs);
+}
+
+void test_random_polygons(CGAL::Random &r, size_t nb_polygons, size_t nb_pts){
+#ifdef BENCH_AND_VERBOSE_FLOAT_SNAP_ROUNDING_2
+  std::cout << "Test random polygons" << std::endl;
+#endif
+  std::vector<Polygon_2> polygons;
+  for(size_t i=0; i!=nb_polygons; ++i){
+    Polygon_2 poly;
+    for(size_t j=0; j!=nb_pts; ++j)
+      poly.push_back(Point_2(random_point(r), random_point(r)));
+    polygons.emplace_back(poly);
+  }
+
+  std::vector<Polygon_2> out;
+  CGAL::compute_snapped_polygons_2(polygons.begin(), polygons.end(), std::back_inserter(out));
+
+  std::vector<Segment_2> segs;
+  for(const Polygon_2 &poly: out)
+    for(size_t i=1; i<poly.size(); ++i)
+      segs.emplace_back(poly[i-1],poly[i]);
+  assert(!CGAL::do_curves_intersect(segs.begin(), segs.end()));
+}
+
+void test_almost_indentical_segments(CGAL::Random &r, size_t nb_segments, Vector_2 source, Vector_2 target){
+  // Difficult test
+#ifdef BENCH_AND_VERBOSE_FLOAT_SNAP_ROUNDING_2
+  std::cout << "Test with almost identical segments from " << source << " to " << target << std::endl;
+#endif
+  std::vector<Segment_2> segs;
+  double eps=std::pow(2,-40);
+  for(size_t i=0; i!=nb_segments; ++i){
+    segs.emplace_back(random_noise_point(r, -eps, eps, source), random_noise_point(r, -eps, eps, target));
+  }
+
+  test(segs);
+}
+
+void test_iterative_square_intersection(CGAL::Random &r, size_t nb_iterations){
+  auto random_rotated_square=[&](){
+    double theta=r.get_double(0, CGAL_PI/2);
+#ifdef BENCH_AND_VERBOSE_FLOAT_SNAP_ROUNDING_2
+    std::cout << "Angle: " << theta << std::endl;
+#endif
+    FT cos_t(std::cos(theta));
+    FT sin_t(std::sin(theta));
+    Polygon_2 P;
+    P.push_back(Point_2( cos_t, sin_t));
+    P.push_back(Point_2(-sin_t, cos_t));
+    P.push_back(Point_2(-cos_t,-sin_t));
+    P.push_back(Point_2( sin_t,-cos_t));
+    return P;
+  };
+
+  Polygon_2 scene=random_rotated_square();
+  Polygon_2 snap_scene;
+  Pwh_vec_2 out_intersection;
+
+  for(size_t i=0; i<nb_iterations; ++i){
+    out_intersection.clear();
+    CGAL::intersection(random_rotated_square(), scene, std::back_inserter(out_intersection));
+    assert(out_intersection.size()==1 && out_intersection[0].number_of_holes()==0);
+#ifdef BENCH_AND_VERBOSE_FLOAT_SNAP_ROUNDING_2
+    CGAL::Real_timer t;
+    t.start();
+#endif
+    compute_snapped_polygon_2(out_intersection[0].outer_boundary(), snap_scene);
+#ifdef BENCH_AND_VERBOSE_FLOAT_SNAP_ROUNDING_2
+    t.stop();
+    std::cout << "Iteration " << i << std::endl;
+    std::cout << "Polygon size: " << out_intersection[0].outer_boundary().size()
+              << " , Snapped polygon size: " << snap_scene.size() << std::endl;
+    std::cout << "is convex: " << snap_scene.is_convex() << std::endl;
+    std::cout << "Running time: " << t.time() << std::endl;
+#endif
+    scene=snap_scene;
+  }
+}
+
+void test_multi_almost_indentical_segments(CGAL::Random &r, size_t nb_segments){
+  for(double x1=-1; x1<=1; ++x1)
+    for(double y1=-1; y1<=1; ++y1)
+      for(double x2=x1; x2<=1; ++x2)
+        for(double y2=y1; y2<=1; ++y2)
+          if(Vector_2(x1,y1)!=Vector_2(x2,y2))
+            test_almost_indentical_segments(r, nb_segments, Vector_2(x1,y1), Vector_2(x2,y2));
+}
+
+void fix_test(){
+#ifdef BENCH_AND_VERBOSE_FLOAT_SNAP_ROUNDING_2
+  std::cout << "Fix tests" << std::endl;
+#endif
+  std::vector< Segment_2 > segs;
+  FT e(std::pow(2, -60));
+
+  // Disjoint segments
+  segs.emplace_back(Point_2(0, 1), Point_2(1, 1));
+  segs.emplace_back(Point_2(1, 2), Point_2(2, 2));
+  test(segs);
+  segs.clear();
+
+  // Degenerate segment
+  segs.emplace_back(Point_2(0, 1), Point_2(2, 1));
+  segs.emplace_back(Point_2(1, 0), Point_2(1, 0));
+  test(segs);
+  segs.clear();
+
+  // Collinear segments
+  segs.emplace_back(Point_2(0, 0), Point_2(4, 0));
+  segs.emplace_back(Point_2(2, 0), Point_2(6, 0));
+  test(segs);
+  segs.clear();
+
+  // Duplicate segment
+  segs.emplace_back(Point_2(0, 0), Point_2(4, 0));
+  segs.emplace_back(Point_2(0, 0), Point_2(4, 0));
+  test(segs);
+  segs.clear();
+
+  // Polyline
+  segs.emplace_back(Point_2(0, 0), Point_2(2, 0));
+  segs.emplace_back(Point_2(2, 0), Point_2(1, 1));
+  segs.emplace_back(Point_2(1, 1), Point_2(3, 1));
+  test(segs);
+  segs.clear();
+
+  // Horizontal intersection when round
+  segs.emplace_back(Point_2(0, 1+e), Point_2(2, 1));
+  segs.emplace_back(Point_2(1, 1), Point_2(1, 0));
+  test(segs);
+  segs.clear();
+
+  // Vertical intersection when round
+  segs.emplace_back(Point_2(0, 0), Point_2(1-e, 0));
+  segs.emplace_back(Point_2(1, 1), Point_2(1, -1));
+  test(segs);
+  segs.clear();
+
+  // A more complex example
+  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));
+  segs.emplace_back(Point_2(0, 2-e/2), Point_2(2, 0));
+  segs.emplace_back(Point_2(0, 2-e/2), Point_2(-2+e, -4));
+  segs.emplace_back(Point_2(-2, 2), Point_2(2, 2));
+  segs.emplace_back(Point_2(7, 7), Point_2(7+e, 7+e));
+  segs.emplace_back(Point_2(5, 7-e), Point_2(9, 7-e));
+  test(segs);
+
+  // Minimal intersection subsets of a random examples
+  segs.emplace_back(Point_2(0.99999999999992173,-0.9999999999991368),Point_2(0.99999999999987266,-6.016984285966173e-13));
+  segs.emplace_back(Point_2(1.000000000000494,-1.0000000000002354),Point_2(0.99999999999950062,8.0391743541535603e-13));
+  segs.emplace_back(Point_2(1.0000000000008489,-0.99999999999951794),Point_2(0.99999999999926925,-2.3642280455149055e-13));
+  test(segs);
+  segs.clear();
+
+  segs.emplace_back(Point_2(1.0000000000008309,9.0061990813884068e-13), Point_2(0.99999999999981259,0.99999999999938394));
+  segs.emplace_back(Point_2(1.0000000000005094,-3.1951552372595695e-13), Point_2(1.0000000000004887,1.0000000000006302));
+  segs.emplace_back(Point_2(1.0000000000006171,-7.3034227828590228e-13), Point_2(1.0000000000004181,0.99999999999913269));
+  segs.emplace_back(Point_2(1.0000000000008491,6.6531536071947998e-14), Point_2(0.9999999999993564,1.0000000000006497));
+  segs.emplace_back(Point_2(1.0000000000006859,3.4155608842536406e-13), Point_2(0.99999999999990252,0.99999999999978451));
+  segs.emplace_back(Point_2(1.0000000000007121,-1.6363168568197086e-13), Point_2(1.0000000000004778,0.99999999999927969));
+  segs.emplace_back(Point_2(1.000000000000753,8.3769047200168284e-13), Point_2(1.0000000000003773,1.0000000000001725));
+  segs.emplace_back(Point_2(1.000000000000574,6.6944274898354285e-13), Point_2(1.000000000000427,0.99999999999998845));
+  segs.emplace_back(Point_2(1.0000000000006153,-2.798953111315494e-13), Point_2(1.0000000000003044,1.0000000000008289));
+  segs.emplace_back(Point_2(1.0000000000008387,5.8214218612051864e-13), Point_2(1.0000000000002918,0.99999999999943212));
+  segs.emplace_back(Point_2(1.0000000000007745,-3.9231414307222458e-13), Point_2(1.000000000000294,1.0000000000003408));
+  test(segs);
+  segs.clear();
+}
+
+void test_float_snap_rounding(){
+#ifdef BENCH_AND_VERBOSE_FLOAT_SNAP_ROUNDING_2
+  std::cout << "Test the other API" << std::endl;
+#endif
+
+  std::vector< Segment_2 > segs;
+  std::vector< Polyline_2 > out;
+  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));
+  segs.emplace_back(Point_2(0, 2-e/2), Point_2(2, 0));
+  segs.emplace_back(Point_2(0, 2-e/2), Point_2(-2+e, -4));
+  segs.emplace_back(Point_2(-2, 2), Point_2(2, 2));
+  segs.emplace_back(Point_2(7, 7), Point_2(7+e, 7+e));
+  segs.emplace_back(Point_2(5, 7-e), Point_2(9, 7-e));
+
+  double_snap_rounding_2(segs.begin(), segs.end(), out);
+}
+
+int main(int argc,char *argv[])
+{
+  CGAL::Random rp;
+  CGAL::Random r(argc==1?rp.get_seed():std::stoi(argv[1]));
+  std::cout << "random seed = " << r.get_seed() << std::endl;
+  std::cout << std::setprecision(17);
+  fix_test();
+  test_float_snap_rounding();
+  test_fully_random(r,1000);
+  test_multi_almost_indentical_segments(r,200);
+  // test_random_polygons(r,200,10);
+  // test_iterative_square_intersection(r,2000);
+  return(0);
+}
diff --git a/Snap_rounding_2/test/Snap_rounding_2/test_float_snap_rounding_2_concept.cpp b/Snap_rounding_2/test/Snap_rounding_2/test_float_snap_rounding_2_concept.cpp
new file mode 100644
index 00000000000..518f8d405e4
--- /dev/null
+++ b/Snap_rounding_2/test/Snap_rounding_2/test_float_snap_rounding_2_concept.cpp
@@ -0,0 +1,42 @@
+// #define DOUBLE_2D_SNAP_VERBOSE
+// #define DOUBLE_2D_SNAP_FULL_VERBOSE
+#define BENCH_AND_VERBOSE_FLOAT_SNAP_ROUNDING_2
+
+#include <CGAL/Exact_rational.h>
+#include <CGAL/Cartesian.h>
+#include <CGAL/Cartesian_converter.h>
+
+#include <CGAL/Float_snap_rounding_traits_2.h>
+#include <CGAL/Float_snap_rounding_2.h>
+
+typedef CGAL::Exact_predicates_exact_constructions_kernel                      Epeck;
+typedef CGAL::Cartesian<CGAL::Exact_rational>                                 Rational_Kernel;
+typedef CGAL::Float_snap_rounding_traits_2<Rational_Kernel, Rational_Kernel>  Rational_Traits;
+
+template<class Traits>
+struct Test{
+  typedef typename Traits::Segment_2 Segment_2;
+  typedef typename Traits::Point_2   Point_2;
+
+  void fix_test(){
+    std::vector<Segment_2> segs;
+    segs.emplace_back(Point_2(1, 1), Point_2(-1, -1));
+    segs.emplace_back(Point_2(0, 0), Point_2(1, -1));
+    segs.emplace_back(Point_2(0, 2), Point_2(2, 0));
+    segs.emplace_back(Point_2(0, 2), Point_2(-2, -4));
+    segs.emplace_back(Point_2(-2, 2), Point_2(2, 2));
+    segs.emplace_back(Point_2(5, 7), Point_2(9, 7));
+    segs.emplace_back(Point_2(1,1), Point_2(3,1));
+    segs.emplace_back(Point_2(1,2), Point_2(3,0));
+    std::vector<Segment_2> out;
+    Traits traits;
+    CGAL::compute_snapped_subcurves_2(segs.begin(), segs.end(), std::back_inserter(out), CGAL::parameters::geom_traits(traits));
+    assert(!CGAL::do_curves_intersect(out.begin(), out.end()));
+  }
+};
+
+int main(/*int argc,char *argv[]*/)
+{
+  Test<Rational_Traits>().fix_test();
+  return(0);
+}
diff --git a/Surface_sweep_2/examples/Surface_sweep_2/plane_sweep.cpp b/Surface_sweep_2/examples/Surface_sweep_2/plane_sweep.cpp
index 28f894cd8bf..1493b1c5f63 100644
--- a/Surface_sweep_2/examples/Surface_sweep_2/plane_sweep.cpp
+++ b/Surface_sweep_2/examples/Surface_sweep_2/plane_sweep.cpp
@@ -6,25 +6,56 @@
 
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 #include <CGAL/Arr_segment_traits_2.h>
+#include <CGAL/Arr_polyline_traits_2.h>
 #include <CGAL/Surface_sweep_2_algorithms.h>
 
 typedef CGAL::Exact_predicates_exact_constructions_kernel       Kernel;
 typedef Kernel::Point_2                                         Point_2;
-typedef CGAL::Arr_segment_traits_2<Kernel>                      Traits_2;
-typedef Traits_2::Curve_2                                       Segment_2;
+typedef Kernel::Segment_2                                         Segment;
+typedef Kernel::Point_2                                         Point;
+using Segment_traits = CGAL::Arr_segment_traits_2<Kernel>;
+using Traits = CGAL::Arr_polyline_traits_2<Segment_traits>;
+using My_polyline = Traits::Curve_2;
 
 int main()
 {
+    Traits traits;
   // Construct the input segments.
-  Segment_2 segments[] = {Segment_2 (Point_2 (1, 5), Point_2 (8, 5)),
-                          Segment_2 (Point_2 (1, 1), Point_2 (8, 8)),
-                          Segment_2 (Point_2 (3, 1), Point_2 (3, 8)),
-                          Segment_2 (Point_2 (8, 5), Point_2 (8, 8))};
+  auto polyline_construct = traits.construct_curve_2_object();
+
+  Point points1[5];
+  points1[0] = Point(0, 0);
+  points1[1] = Point(2, 4);
+  points1[2] = Point(3, 0);
+  points1[3] = Point(4, 4);
+  points1[4] = Point(6, 0);
+  auto pi1 = polyline_construct(&points1[0], &points1[5]);
+
+  std::list<Point> points2;
+  points2.push_back(Point(1, 3));
+  points2.push_back(Point(0, 2));
+  points2.push_back(Point(1, 0));
+  points2.push_back(Point(2, 1));
+  points2.push_back(Point(3, 0));
+  points2.push_back(Point(4, 1));
+  points2.push_back(Point(5, 0));
+  points2.push_back(Point(6, 2));
+  points2.push_back(Point(5, 3));
+  points2.push_back(Point(4, 2));
+  auto pi2 = polyline_construct(points2.begin(), points2.end());
+
+  std::vector<Segment> segs;
+  segs.push_back(Segment(Point(0, 2), Point(1, 2)));
+  segs.push_back(Segment(Point(1, 2), Point(3, 6)));
+  segs.push_back(Segment(Point(3, 6), Point(5, 2)));
+  auto pi3 = polyline_construct(segs.begin(), segs.end());
+
+  My_polyline vec[3]={pi1,  pi2, pi3};
 
   // Compute all intersection points.
   std::list<Point_2> pts;
 
-  CGAL::compute_intersection_points(segments, segments + 4,
+  CGAL::compute_intersection_points(vec, vec + 3,
                                     std::back_inserter(pts));
 
   // Print the result.
@@ -33,14 +64,14 @@ int main()
             std::ostream_iterator<Point_2>(std::cout, "\n"));
 
   // Compute the non-intersecting sub-segments induced by the input segments.
-  std::list<Segment_2> sub_segs;
+  std::list<My_polyline> sub_segs;
 
-  CGAL::compute_subcurves(segments, segments + 4, std::back_inserter(sub_segs));
+  CGAL::compute_subcurves(vec, vec + 3, std::back_inserter(sub_segs));
 
   std::cout << "Found " << sub_segs.size()
             << " interior-disjoint sub-segments." << std::endl;
 
-  assert(CGAL::do_curves_intersect (segments, segments + 4));
+//  assert(CGAL::do_curves_intersect (segments, segments + 4));
 
   return 0;
 }
diff --git a/Surface_sweep_2/include/CGAL/Surface_sweep_2.h b/Surface_sweep_2/include/CGAL/Surface_sweep_2.h
index 2c007c5955a..93d46ecc8a9 100644
--- a/Surface_sweep_2/include/CGAL/Surface_sweep_2.h
+++ b/Surface_sweep_2/include/CGAL/Surface_sweep_2.h
@@ -198,7 +198,7 @@ protected:
 
   /*! Compute intersections between the two given curves.
    * If the two curves intersect, create a new event (or use the event that
-   * already exits in the intersection point) and insert the curves to the
+   * already exists in the intersection point) and insert the curves to the
    * event.
    * \param curve1 The first curve.
    * \param curve2 The second curve.