diff --git a/Installation/CHANGES.md b/Installation/CHANGES.md
index d802f740b9b..09bde261d24 100644
--- a/Installation/CHANGES.md
+++ b/Installation/CHANGES.md
@@ -14,7 +14,7 @@ Release date: October 2023
 - **Breaking change**: The usage of `boost::shared_ptr` has been replaced by `std::shared_ptr`. Packages affected are 2D Straight Line Skeleton and Shape Detection.
 - **Breaking change**: The usage of `boost::optional` has been replaced by `std::optional`. Packages affected are 2D Straight Line Skeleton, 3D Fast Intersection and Distance Computation (AABB Tree), and the Kernel intersection.
 - **Breaking change**: The usage of `boost::variant` has been replaced by `std::variant`. Packages affected are 2D Arrangements, and the Kernel intersection.
-- **Breaking chahge**: The file CMake file `UseCGAL.cmake` has been removed from CGAL. Usages of the CMake variables `${CGAL_USE_FILE}` and `${CGAL_LIBRARIES}` must be replaced by a link to the imported target `CGAL::CGAL`, for example: `target_link_library(the_target PRIVATE CGAL::CGAL)`.
+- **Breaking change**: The file CMake file `UseCGAL.cmake` has been removed from CGAL. Usages of the CMake variables `${CGAL_USE_FILE}` and `${CGAL_LIBRARIES}` must be replaced by a link to the imported target `CGAL::CGAL`, for example: `target_link_library(the_target PRIVATE CGAL::CGAL)`.
 
 
 #### 2D Arrangements
@@ -45,6 +45,10 @@ Release date: October 2023
 -   Added the option to use a variable sizing field for `CGAL::Polygon_mesh_processing::isotropic_remeshing()`,
     and a sizing function based on a measure of local curvature for adaptive remeshing.
 -   Added the function `CGAL::Polygon_mesh_processing::refine_mesh_at_isolevel()` that refines a polygon mesh along an isocurve.
+-   Added the function
+    `CGAL::Polygon_mesh_processing::autorefine_triangle_soup()` that refines a soup of triangles so that no pair of triangles intersects
+     in their interiors. Also added, the function `autorefine()` operating directly on a triangle mesh and updating it
+     using the aforementioned function on a triangle soup.
 
 ### [2D Arrangements](https://doc.cgal.org/6.0/Manual/packages.html#PkgArrangementOnSurface2)
 -   Fixed a bug in the zone construction code applied to arrangements of geodesic arcs on a sphere,
diff --git a/Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Triangle_3_intersection.h b/Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Triangle_3_intersection.h
index 42182e01636..269c86c3e60 100644
--- a/Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Triangle_3_intersection.h
+++ b/Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Triangle_3_intersection.h
@@ -14,6 +14,8 @@
 #ifndef CGAL_INTERNAL_INTERSECTIONS_TRIANGLE_3_TRIANGLE_3_INTERSECTION_H
 #define CGAL_INTERNAL_INTERSECTIONS_TRIANGLE_3_TRIANGLE_3_INTERSECTION_H
 
+//#define CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+
 #include <CGAL/Intersection_traits_3.h>
 #include <CGAL/Intersections_3/internal/Line_3_Triangle_3_intersection.h>
 #include <CGAL/Intersections_3/internal/Line_3_Line_3_intersection.h>
@@ -21,6 +23,9 @@
 
 #include <CGAL/kernel_assertions.h>
 
+#include <boost/next_prior.hpp>
+#include <boost/container/flat_set.hpp>
+
 #include <list>
 #include <map>
 #include <vector>
@@ -30,53 +35,338 @@ namespace Intersections {
 namespace internal{
 
 template <class Kernel>
-void intersection_coplanar_triangles_cutoff(const typename Kernel::Point_3& p,
-                                            const typename Kernel::Point_3& q,
-                                            const typename Kernel::Point_3& r,
-                                            const Kernel& k,
-                                            std::list<typename Kernel::Point_3>& inter_pts)
+struct Point_on_triangle
 {
-  typedef typename std::list<typename Kernel::Point_3>::iterator Iterator;
+  // triangle points are not stored in this class but are expected
+  // to always be passed in the same order. For a triangle pqr,
+  // edge 0 is pq, edge 1 qr and edge 2 rp. Point 0 is p, 1 is q and 2 is r.
+  //
+  // (id, -1) point on t1
+  // (-1, id) point on t2
+  // (id1, id2) intersection of edges
+  std::pair<int, int> t1_t2_ids;
+  boost::container::flat_set<int> extra_t1; // store other ids of edges containing the point
+  typename Kernel::FT alpha; //
+
+//////
+
+  static
+  inline
+  const typename Kernel::Point_3&
+  point_from_id(const typename Kernel::Point_3& p,
+                const typename Kernel::Point_3& q,
+                const typename Kernel::Point_3& r,
+                int id)
+  {
+    switch(id)
+    {
+      case 0:
+        return p;
+      case 1:
+        return q;
+      default:
+        return r;
+    }
+  }
+
+  Point_on_triangle(int i1, int i2=-1, typename Kernel::FT alpha = 0.) // TODO add global zero()?
+    : t1_t2_ids(i1,i2)
+    , alpha(alpha)
+  {}
+
+  // orientation of the current point wrt to edge id1 (p1q1)
+  Orientation
+  orientation (const typename Kernel::Point_3& p1, // source of edge edge_id1
+               const typename Kernel::Point_3& q1, // target of edge edge_id1
+               const typename Kernel::Point_3& r1,
+               int edge_id1,
+               const typename Kernel::Point_3& p2,
+               const typename Kernel::Point_3& q2,
+               const typename Kernel::Point_3& r2,
+               const Kernel& k) const
+  {
+    if (t1_t2_ids.first!=-1)
+    {
+      if (t1_t2_ids.second==-1)
+        return (edge_id1==t1_t2_ids.first || (edge_id1+1)%3==t1_t2_ids.first) ? ZERO:POSITIVE; // it is a point on t1
+      // this is an intersection point
+
+      if (t1_t2_ids.first==edge_id1)
+        return ZERO;
+      if (t1_t2_ids.first==(edge_id1+1)%3)
+      {
+        if (alpha==0) return ZERO;
+        return alpha>=0 ? POSITIVE:NEGATIVE;
+      }
+      CGAL_assertion((t1_t2_ids.first+1)%3==edge_id1);
+      if (alpha==1) return ZERO;
+      return alpha<=1?POSITIVE:NEGATIVE;
+    }
+    else
+    {
+      //this is an input point of t2
+      typename Kernel::Coplanar_orientation_3 orient = k.coplanar_orientation_3_object();
+      const typename Kernel::Point_3& query = point_from_id(p2,q2,r2,t1_t2_ids.second);
+      return orient(p1,q1,r1,query);
+    }
+  }
+
+  int id1() const { return t1_t2_ids.first; }
+  int id2() const { return t1_t2_ids.second; }
+
+  // construct the intersection point from the info stored
+  typename Kernel::Point_3
+  point(const typename Kernel::Point_3& p1,
+        const typename Kernel::Point_3& q1,
+        const typename Kernel::Point_3& r1,
+        const typename Kernel::Point_3& p2,
+        const typename Kernel::Point_3& q2,
+        const typename Kernel::Point_3& r2,
+        const Kernel& k) const
+  {
+    if (t1_t2_ids.first==-1)
+      return point_from_id(p2,q2,r2,t1_t2_ids.second);
+    if (t1_t2_ids.second==-1)
+      return point_from_id(p1,q1,r1,t1_t2_ids.first);
+
+    return k.construct_barycenter_3_object()(point_from_id(p1,q1,r1,(t1_t2_ids.first+1)%3), alpha, point_from_id(p1,q1,r1,t1_t2_ids.first)) ;
+  }
+};
+
+// the intersection of two triangles is computed by interatively intersection t2
+// with halfspaces defined by edges of t1. The following function is called
+// for each each on t1 on edge of the current intersection.
+// pq is such an edge and p1q1 from t1 defines the halfspace intersection
+// we are currently interseted in. We return the intersection point of
+// pq with p1q1
+template <class Kernel>
+Point_on_triangle<Kernel>
+intersection(const Point_on_triangle<Kernel>& p,
+             const Point_on_triangle<Kernel>& q,
+             int edge_id_t1,
+             const typename Kernel::Point_3& p1,
+             const typename Kernel::Point_3& q1,
+//             const typename Kernel::Point_3& r1,
+             const typename Kernel::Point_3& p2,
+             const typename Kernel::Point_3& q2,
+             const typename Kernel::Point_3& r2,
+             const Kernel& k)
+{
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+  std::cout << "    calling intersection: ";
+  std::cout << " (" << p.id1() << "," << p.id2() << ",[" << p.alpha << "]) -";
+  std::cout << " (" << q.id1() << "," << q.id2() << ",[" << q.alpha << "]) || e" << edge_id_t1;
+#endif
+  typename Kernel::Compute_alpha_for_coplanar_triangle_intersection_3 compute_alpha
+    = k.compute_alpha_for_coplanar_triangle_intersection_3_object();
+  typedef Point_on_triangle<Kernel> Pot;
+  switch(p.id1())
+  {
+    case -1:
+    {
+      switch(q.id1())
+      {
+        case -1: // A: (-1, ip2) - (-1, iq2)
+        {
+          CGAL_assertion((p.id2()+1)%3 == q.id2() || (q.id2()+1)%3 == p.id2());
+//          CGAL_assertion(p.extra_t1.empty() && q.extra_t1.empty());  // TMP to see if it's worth implementing special case
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+          std::cout << " -- case 1\n";
+#endif
+          typename Kernel::FT alpha = compute_alpha(p1, q1,
+                                                    Pot::point_from_id(p2, q2, r2, p.id2()),
+                                                    Pot::point_from_id(p2, q2, r2, q.id2()));
+          int id2 = (p.id2()+1)%3 == q.id2() ? p.id2() : q.id2();
+          return  Point_on_triangle<Kernel>(edge_id_t1, id2, alpha); // intersection with an original edge of t2
+        }
+        default:
+          if (q.id2()!=-1) // B: (-1, ip2) - (iq1, iq2)
+          {
+            if (p.id2() == q.id2() || p.id2() == (q.id2()+1)%3)
+            {
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+              std::cout << " -- case 2\n";
+#endif
+              // points are on the same edge of t2 --> we shorten an already cut edge
+              typename Kernel::FT alpha = compute_alpha(p1, q1,
+                                                        Pot::point_from_id(p2, q2, r2,  q.id2()),
+                                                        Pot::point_from_id(p2, q2, r2, (q.id2()+1)%3));
+
+              return  Point_on_triangle<Kernel>(edge_id_t1, q.id2(), alpha);
+            }
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+            std::cout << " -- case 3\n";
+#endif
+            // point of t1: look for an edge of t1 containing both points
+            CGAL_assertion( p.extra_t1.count(q.id1())!=0 || p.extra_t1.count(3-q.id1()-edge_id_t1)!=0 );
+            int eid1 = p.extra_t1.count(q.id1())!=0 ? q.id1() : 3-q.id1()-edge_id_t1;
+            return  Point_on_triangle<Kernel>((eid1+1)%3==edge_id_t1?edge_id_t1:(edge_id_t1+1)%3, -1); // vertex of t1
+          }
+          // C: (-1, ip2) - (iq1, -1)
+          //vertex of t1, special case t1 edge passed thru a vertex of t2
+          CGAL_assertion(edge_id_t1 == 2);
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+          std::cout << " -- case 4\n";
+#endif
+          CGAL_assertion(q.id1()==1);
+          CGAL_assertion(!p.extra_t1.empty());
+          return  Point_on_triangle<Kernel>(p.extra_t1.count(0)==1?0:2,-1);
+      }
+    }
+    default:
+    {
+      switch(p.id2())
+      {
+        case -1:
+        {
+          switch(q.id1())
+          {
+            case -1: // G: (ip1, -1) - (-1, iq2)
+              //vertex of t1, special case t1 edge passed thru a vertex of t2
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+              std::cout << " -- case 5\n";
+#endif
+              CGAL_assertion(edge_id_t1 == 2);
+              CGAL_assertion(p.id1()==1);
+              CGAL_assertion(!q.extra_t1.empty());
+              return  Point_on_triangle<Kernel>(q.extra_t1.count(0)==1?0:2,-1);
+            default:
+            {
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+              std::cout << " -- case 6\n";
+#endif
+              CGAL_assertion(q.id2()!=-1); // I: (ip1, -1) - (iq2, -1)
+              //H: (ip1,-1), (iq1, iq2)
+              CGAL_assertion(edge_id_t1==2);
+              // p and q are on the same edge of t1
+              CGAL_assertion(p.id1()==q.id1() || p.id1()==(q.id1()+1)%3);
+              return  Point_on_triangle<Kernel>((q.id1()+1)%3==edge_id_t1?edge_id_t1:(edge_id_t1+1)%3 , -1);
+            }
+          }
+        }
+        default:
+        {
+          switch(q.id1())
+          {
+            case -1: // D: (ip1, ip2) - (-1, iq2)
+            {
+              if (q.id2() == p.id2() || q.id2() == (p.id2()+1)%3)
+              {
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+                std::cout << " -- case 7\n";
+#endif
+                // points are on the same edge of t2 --> we shorten an already cut edge
+                typename Kernel::FT alpha = compute_alpha(p1, q1,
+                                                          Pot::point_from_id(p2, q2, r2,  p.id2()),
+                                                          Pot::point_from_id(p2, q2, r2, (p.id2()+1)%3));
+
+                return  Point_on_triangle<Kernel>(edge_id_t1, p.id2(), alpha);
+              }
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+              std::cout << " -- case 8\n";
+#endif
+              // point of t1
+              //std::cout << "q.extra_t1: "; for(int qet1 : q.extra_t1) std::cout << " " << qet1; std::cout << "\n";
+              CGAL_assertion( q.extra_t1.count(p.id1())!=0 || q.extra_t1.count(3-p.id1()-edge_id_t1)!=0 );
+              int eid1 = q.extra_t1.count(p.id1())!=0 ? p.id1() : 3-p.id1()-edge_id_t1;
+              return  Point_on_triangle<Kernel>((eid1+1)%3==edge_id_t1?edge_id_t1:(edge_id_t1+1)%3, -1); // vertex of t1
+            }
+            default:
+            {
+              switch(q.id2())
+              {
+                case -1: // F: (ip1, ip2) - (iq1, -1)
+                {
+                  // p and q are on the same edge of t1
+                  CGAL_assertion(q.id1()==p.id1() ||  q.id1()==(p.id1()+1)%3);
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+                  std::cout << " -- case 9\n";
+#endif
+                  return  Point_on_triangle<Kernel>((p.id1()+1)%3==edge_id_t1?edge_id_t1:(edge_id_t1+1)%3 , -1);
+                }
+                default: // E: (ip1, ip2) - (iq1, iq2)
+                {
+                  if (p.id2()==q.id2())
+                  {
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+                    std::cout << " -- case 10\n";
+#endif
+                    typename Kernel::FT alpha = compute_alpha(p1, q1,
+                                                              Pot::point_from_id(p2, q2, r2,  q.id2()),
+                                                              Pot::point_from_id(p2, q2, r2, (q.id2()+1)%3));
+                    return  Point_on_triangle<Kernel>(edge_id_t1, q.id2(), alpha);
+                  }
+                  // we are intersecting an edge of t1
+                  CGAL_assertion(p.id1()==q.id1() || edge_id_t1==2);
+                  int eid1 = p.id1()==q.id1() ? p.id1() : 1;
+                  return  Point_on_triangle<Kernel>((eid1+1)%3==edge_id_t1?edge_id_t1:(edge_id_t1+1)%3, -1); // vertex of t1
+                }
+              }
+            }
+          }
+        }
+      }
+    }
+  }
+}
+
+template <class Kernel>
+void intersection_coplanar_triangles_cutoff(const typename Kernel::Point_3& p1,
+                                            const typename Kernel::Point_3& q1,
+                                            const typename Kernel::Point_3& r1,
+                                            int edge_id,
+                                            const typename Kernel::Point_3& p2,
+                                            const typename Kernel::Point_3& q2,
+                                            const typename Kernel::Point_3& r2,
+                                            const Kernel& k,
+                                            std::list<Point_on_triangle<Kernel>>& inter_pts)
+{
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+  std::cout << "  cutoff using e" << edge_id << ": "
+            << to_double(p1.x()) << " " << to_double(p1.y()) << " " << to_double(p1.z()) << " "
+            << to_double(q1.x()) << " " << to_double(q1.y()) << " " << to_double(q1.z()) << "\n";
+#endif
+  typedef typename std::list<Point_on_triangle<Kernel>>::iterator Iterator;
 
   if(inter_pts.empty())
     return;
 
-  typename Kernel::Coplanar_orientation_3 orient = k.coplanar_orientation_3_object();
-  typename Kernel::Construct_line_3 line = k.construct_line_3_object();
+  //orient(p1,q1,r1,r1) is POSITIVE
+  std::map<const Point_on_triangle<Kernel>*,Orientation> orientations; // TODO skip map
+  for (Point_on_triangle<Kernel>& pot : inter_pts)
+  {
+    orientations[ &pot ]=pot.orientation(p1,q1,r1,edge_id,p2,q2,r2,k);
+    if (pot.id1()==-1 && orientations[ &pot ]==COLLINEAR)
+      pot.extra_t1.insert(edge_id);
+  }
 
-  //orient(p,q,r,r) is POSITIVE
-  std::map<const typename Kernel::Point_3*,Orientation> orientations;
-  for (Iterator it=inter_pts.begin();it!=inter_pts.end();++it)
-    orientations[ &(*it) ]=orient(p,q,r,*it);
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+  std::cout << "    Orientations:";
+  for (const Point_on_triangle<Kernel>& pot : inter_pts)
+    std::cout << " " << orientations[ &pot ];
+  std::cout << "\n";
+#endif
+  CGAL_kernel_assertion_code(int pt_added = 0);
 
-  CGAL_kernel_assertion_code(int pt_added = 0;)
+  Iterator prev = std::prev(inter_pts.end());
 
-  const typename Kernel::Point_3* prev = &(*std::prev(inter_pts.end()));
   Iterator stop = inter_pts.size() > 2 ? inter_pts.end() : std::prev(inter_pts.end());
   for(Iterator it=inter_pts.begin(); it!=stop; ++it)
   {
-    const typename Kernel::Point_3& curr = *it;
-    Orientation or_prev = orientations[prev],
-                or_curr = orientations[&curr];
+    Orientation or_prev = orientations[&(*prev)],
+                or_curr = orientations[&(*it)];
 
     if((or_prev == POSITIVE && or_curr == NEGATIVE) ||
        (or_prev == NEGATIVE && or_curr == POSITIVE))
     {
-      typename Intersection_traits<Kernel, typename Kernel::Line_3, typename Kernel::Line_3>::result_type
-        obj = intersection(line(p,q), line(*prev,curr), k);
+      Point_on_triangle<Kernel> new_pt = intersection(*prev, *it, edge_id, p1, q1, p2, q2, r2, k);
 
-      // assert "not empty"
-      CGAL_kernel_assertion(bool(obj));
-
-      const typename Kernel::Point_3* inter = intersect_get<typename Kernel::Point_3>(obj);
-      CGAL_kernel_assertion(inter != nullptr);
-
-      prev = &(*inter_pts.insert(it,*inter));
-      orientations[prev] = COLLINEAR;
-      CGAL_kernel_assertion_code(++pt_added;)
+      prev = inter_pts.insert(it,new_pt);
+      orientations[&(*prev)] = COLLINEAR;
+      CGAL_kernel_assertion_code(++pt_added);
     }
 
-    prev = &(*it);
+    prev = it;
   }
 
   CGAL_kernel_assertion(pt_added<3);
@@ -96,35 +386,77 @@ intersection_coplanar_triangles(const typename K::Triangle_3& t1,
                                 const typename K::Triangle_3& t2,
                                 const K& k)
 {
-  const typename K::Point_3& p = t1.vertex(0),
-                             q = t1.vertex(1),
-                             r = t1.vertex(2);
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+  auto to_string = [](const typename K::Triangle_3& t)
+  {
+    std::stringstream sstr;
+    sstr << "4 "
+         << to_double(t[0].x()) << " "  << to_double(t[0].y()) << " "  << to_double(t[0].z()) << " "
+         << to_double(t[1].x()) << " "  << to_double(t[1].y()) << " "  << to_double(t[1].z()) << " "
+         << to_double(t[2].x()) << " "  << to_double(t[2].y()) << " "  << to_double(t[2].z()) << " "
+         << to_double(t[0].x()) << " "  << to_double(t[0].y()) << " "  << to_double(t[0].z()) << "\n";
+    return sstr.str();
+  };
 
-  std::list<typename K::Point_3> inter_pts;
-  inter_pts.push_back(t2.vertex(0));
-  inter_pts.push_back(t2.vertex(1));
-  inter_pts.push_back(t2.vertex(2));
+  std::cout << "intersection_coplanar_triangles\n";
+  std::ofstream("/tmp/t1.polylines.txt") << std::setprecision(17) << to_string(t1) << "\n";
+  std::ofstream("/tmp/t2.polylines.txt") << std::setprecision(17) << to_string(t2) << "\n";
+#endif
+  const typename K::Point_3& p1 = t1.vertex(0),
+                             q1 = t1.vertex(1),
+                             r1 = t1.vertex(2);
 
+  const typename K::Point_3& p2 = t2.vertex(0),
+                             q2 = t2.vertex(1),
+                             r2 = t2.vertex(2);
+
+  std::list<Point_on_triangle<K>> inter_pts;
+  inter_pts.push_back(Point_on_triangle<K>(-1,0));
+  inter_pts.push_back(Point_on_triangle<K>(-1,1));
+  inter_pts.push_back(Point_on_triangle<K>(-1,2));
+
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+  auto print_points = [&]()
+  {
+    for(auto p : inter_pts) std::cout << "  (" << p.id1() << "," << p.id2() << ",[" << p.alpha << "]) "; std::cout <<"\n";
+  };
+  std::cout << "  ipts size: " << inter_pts.size() << "\n";
+  print_points();
+#endif
   //intersect t2 with the three half planes which intersection defines t1
-  intersection_coplanar_triangles_cutoff(p,q,r,k,inter_pts); //line pq
-  intersection_coplanar_triangles_cutoff(q,r,p,k,inter_pts); //line qr
-  intersection_coplanar_triangles_cutoff(r,p,q,k,inter_pts); //line rp
+  intersection_coplanar_triangles_cutoff(p1,q1,r1,0,p2,q2,r2,k,inter_pts); //line pq
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+  std::cout << "  ipts size: " << inter_pts.size() << "\n";
+  print_points();
+#endif
+  intersection_coplanar_triangles_cutoff(q1,r1,p1,1,p2,q2,r2,k,inter_pts); //line qr
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+  std::cout << "  ipts size: " << inter_pts.size() << "\n";
+  print_points();
+#endif
+  intersection_coplanar_triangles_cutoff(r1,p1,q1,2,p2,q2,r2,k,inter_pts); //line rp
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+  std::cout << "  ipts size: " << inter_pts.size() << "\n";
+  print_points();
+#endif
 
+  auto point = [&](const Point_on_triangle<K>& pot){ return pot.point(p1,q1,r1,p2,q2,r2,k); };
   switch(inter_pts.size())
   {
     case 0:
       return intersection_return<typename K::Intersect_3, typename K::Triangle_3, typename K::Triangle_3>();
     case 1:
-      return intersection_return<typename K::Intersect_3, typename K::Triangle_3, typename K::Triangle_3>(*inter_pts.begin());
+      return intersection_return<typename K::Intersect_3, typename K::Triangle_3, typename K::Triangle_3>(point(*inter_pts.begin()));
     case 2:
       return intersection_return<typename K::Intersect_3, typename K::Triangle_3, typename K::Triangle_3>(
-            k.construct_segment_3_object()(*inter_pts.begin(), *std::next(inter_pts.begin())) );
+            k.construct_segment_3_object()(point(*inter_pts.begin()), point(*std::next(inter_pts.begin()))) );
     case 3:
       return intersection_return<typename K::Intersect_3, typename K::Triangle_3, typename K::Triangle_3>(
-            k.construct_triangle_3_object()(*inter_pts.begin(), *std::next(inter_pts.begin()), *std::prev(inter_pts.end())) );
+            k.construct_triangle_3_object()(point(*inter_pts.begin()), point(*std::next(inter_pts.begin())), point(*std::prev(inter_pts.end()))) );
     default:
       return intersection_return<typename K::Intersect_3, typename K::Triangle_3, typename K::Triangle_3>(
-            std::vector<typename K::Point_3>(inter_pts.begin(),inter_pts.end()));
+            std::vector<typename K::Point_3>(boost::make_transform_iterator(inter_pts.begin(), point),
+                                             boost::make_transform_iterator(inter_pts.end(), point)));
   }
 }
 
diff --git a/Intersections_3/test/Intersections_3/triangle_3_triangle_3_intersection.cpp b/Intersections_3/test/Intersections_3/triangle_3_triangle_3_intersection.cpp
index 0098bcd0142..f0875a55ee8 100644
--- a/Intersections_3/test/Intersections_3/triangle_3_triangle_3_intersection.cpp
+++ b/Intersections_3/test/Intersections_3/triangle_3_triangle_3_intersection.cpp
@@ -96,6 +96,13 @@ void test_coplanar_triangles(){
   assert(CGAL::object_cast<Triangle>(&obj)!=nullptr);
   obj=CGAL::intersection(t2,t1);
   assert(CGAL::object_cast<Triangle>(&obj)!=nullptr);
+  // TK10 case C'
+  t1=Triangle(Point(88.7921, 89.0007, 1.25), Point(88.1912, 88.3997, 1.25), Point(89.8224, 90.031, 1.25));
+  t2=Triangle(Point(88.0497, 88.2583, 1.25), Point(82.9292, 81.8747, 1.25), Point(91.1726, 91.3812, 1.25));
+  obj=CGAL::intersection(t1,t2);
+  assert(CGAL::object_cast<Triangle>(&obj)!=nullptr);
+  obj=CGAL::intersection(t2,t1);
+  assert(CGAL::object_cast<Triangle>(&obj)!=nullptr);
 //Intersection is a point
   //edges  are collinear, one vertex in common
   t1=Triangle( Point(0,0,0),Point(0,1,0),Point(1,0,0) );
@@ -153,6 +160,13 @@ void test_coplanar_triangles(){
   assert(CGAL::object_cast<Segment>(&obj)!=nullptr);
   obj=CGAL::intersection(t2,t1);
   assert(CGAL::object_cast<Segment>(&obj)!=nullptr);
+  // TK10 case D
+  t1=Triangle(Point(-34.893700000000003, -16.0351, 3.1334899999999998e-12), Point(-34.893700000000003, -18.5351, 3.1334899999999998e-12), Point(-42.393700000000003, -16.0351, 3.1334899999999998e-12));
+  t2=Triangle(Point(-34.893700000000003, -32.0351, 3.1334899999999998e-12), Point(-34.893700000000003, -9.7851400000000002, 3.1334899999999998e-12), Point(-31.643699999999999, -17.201799999999999, 3.1334899999999998e-12));
+  obj=CGAL::intersection(t1,t2);
+  assert(CGAL::object_cast<Segment>(&obj)!=nullptr);
+  obj=CGAL::intersection(t2,t1);
+  assert(CGAL::object_cast<Segment>(&obj)!=nullptr);
 //Intersection is a polygon
   //David's star
   t1=Triangle( Point(0,0,0),Point(1,0,0),Point(0.5,1.5,0) );
@@ -181,6 +195,51 @@ void test_coplanar_triangles(){
   obj=CGAL::intersection(t2,t1);
   assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
   assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
+  // TK10 case A
+  t1=Triangle(Point(3.74861, 12.4822, 14.0112), Point(5.40582, 12.4822, 15.6895), Point(5.37748, 12.4822, 15.7206));
+  t2=Triangle(Point(5.49972, 12.4822, 13.491), Point(5.27627, 12.4822, 15.8106), Point(5.32119, 12.4822, 15.8126));
+  obj=CGAL::intersection(t1,t2);
+  assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
+  assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
+  obj=CGAL::intersection(t2,t1);
+  assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
+  assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
+  // TK10 case C
+  t1=Triangle(Point(5, -94.6659, 3.85175), Point(5, -94.5682, 3.08638), Point(5, -94.8182, 3.08638));
+  t2=Triangle(Point(5, -94.4317, 3.76399), Point(5, -97.6182, 3.08638), Point(5, -94.5659, 2.99682));
+  obj=CGAL::intersection(t1,t2);
+  assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
+  assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
+  obj=CGAL::intersection(t2,t1);
+  assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
+  assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
+  // TK10 case E
+  t1=Triangle(Point(-955.858, -45.032, -0.016), Point(-955.856, -45.032, -0.004), Point(-955.856, -45.032, -0.002));
+  t2=Triangle(Point(-955.856, -45.032, 0.006), Point(-955.854, -45.032, -0.002), Point(-955.876, -45.032, -0.034));
+  obj=CGAL::intersection(t1,t2);
+  assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
+  assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
+  obj=CGAL::intersection(t2,t1);
+  assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
+  assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
+  // TK10 case F
+  t1=Triangle(Point(141.172, 20.576, 155.764), Point(141.172, 20.588, 155.766), Point(141.172, 20.59, 155.766));
+  t2=Triangle(Point(141.172, 20.602, 155.768), Point(141.172, 20.594, 155.766), Point(141.172, 20.574, 155.764));
+  obj=CGAL::intersection(t1,t2);
+  assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
+  assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
+  obj=CGAL::intersection(t2,t1);
+  assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
+  assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
+  // TK10 case D
+  t1=Triangle(Point(152.864, 126.324, 0.950001), Point(152.77, 126.483, 0.950001), Point(153.072, 125.973, 0.950001));
+  t2=Triangle(Point(153.322, 125.551, 0.950001), Point(152.218, 127.415, 0.950001), Point(153.66, 124.768, 0.950001));
+  obj=CGAL::intersection(t1,t2);
+  assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
+  assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
+  obj=CGAL::intersection(t2,t1);
+  assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
+  assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
 //Intersection is empty
   t1=Triangle( Point(0,0,0),Point(0,1,0),Point(1,0,0) );
   t2=Triangle( Point(-0.1,-0.1,0),Point(-0.1,-0.9,0),Point(-1,-0.1,0) );
diff --git a/Kernel_23/include/CGAL/Kernel/function_objects.h b/Kernel_23/include/CGAL/Kernel/function_objects.h
index 29a72d9c985..084a6578913 100644
--- a/Kernel_23/include/CGAL/Kernel/function_objects.h
+++ b/Kernel_23/include/CGAL/Kernel/function_objects.h
@@ -2224,6 +2224,107 @@ namespace CommonKernelFunctors {
     }
   };
 
+  template <typename K>
+  class Construct_planes_intersection_point_3
+  {
+    typedef typename K::Plane_3 Plane;
+    typedef typename K::Point_3 Point;
+    typename K::Construct_plane_3 construct_plane;
+  public:
+    typedef Point result_type;
+
+    Point
+    operator()(const Point& p1, const Point& q1, const Point& r1,
+               const Point& p2, const Point& q2, const Point& r2,
+               const Point& p3, const Point& q3, const Point& r3) const
+    {
+      Plane plane1 = construct_plane(p1, q1, r1);
+      Plane plane2 = construct_plane(p2, q2, r2);
+      Plane plane3 = construct_plane(p3, q3, r3);
+
+      const auto res = typename K::Intersect_3()(plane1, plane2, plane3);
+      CGAL_assertion(res!=std::nullopt);
+      const Point* e_pt = std::get_if<Point>(&(*res));
+      CGAL_assertion(e_pt!=nullptr);
+      return *e_pt;
+    }
+
+    Point
+    operator()(const Plane& plane1, const Plane& plane2, const Plane& plane3) const
+    {
+      const auto res = typename K::Intersect_3()(plane1, plane2, plane3);
+      CGAL_assertion(res!=std::nullopt);
+      const Point* e_pt = std::get_if<Point>(&(*res));
+      CGAL_assertion(e_pt!=nullptr);
+      return *e_pt;
+    }
+  };
+
+  template <typename K>
+  class Construct_coplanar_segments_intersection_point_3
+  {
+    typedef typename K::Segment_3 Segment;
+    typedef typename K::Point_3 Point;
+    typename K::Construct_segment_3 construct_segment;
+  public:
+    typedef Point result_type;
+
+    Point
+    operator()(const Point& p1, const Point& q1,
+               const Point& p2, const Point& q2) const
+    {
+      Segment s1 = construct_segment(p1, q1);
+      Segment s2 = construct_segment(p2, q2);
+
+      const auto res = typename K::Intersect_3()(s1, s2);
+      CGAL_assertion(res!=std::nullopt);
+      const Point* e_pt = std::get_if<Point>(&(*res));
+      CGAL_assertion(e_pt!=nullptr);
+      return *e_pt;
+    }
+
+    Point
+    operator()(const Segment& s1, const Segment& s2) const
+    {
+      const auto res = typename K::Intersect_3()(s1, s2);
+      CGAL_assertion(res!=std::nullopt);
+      const Point* e_pt = std::get_if<Point>(&(*res));
+      CGAL_assertion(e_pt!=nullptr);
+      return *e_pt;
+    }
+  };
+
+  template <typename K>
+  class Compute_alpha_for_coplanar_triangle_intersection_3
+  {
+    typedef typename K::Point_3 Point_3;
+    typedef typename K::Vector_3 Vector_3;
+  public:
+    typedef typename K::FT  result_type;
+    result_type
+    operator()(const Point_3& p1, const Point_3& p2,       // segment 1
+               const Point_3& p3, const Point_3& p4) const // segment 2
+    {
+      typename K::Construct_vector_3 vector = K().construct_vector_3_object();
+      typename K::Construct_cross_product_vector_3 cross_product =
+        K().construct_cross_product_vector_3_object();
+
+      const Vector_3 v1 = vector(p1, p2);
+      const Vector_3 v2 = vector(p3, p4);
+
+      CGAL_assertion(K().coplanar_3_object()(p1,p2,p3,p4));
+
+      const Vector_3 v3 = vector(p1, p3);
+      const Vector_3 v3v2 = cross_product(v3,v2);
+      const Vector_3 v1v2 = cross_product(v1,v2);
+      const typename K::FT sl = K().compute_squared_length_3_object()(v1v2);
+      CGAL_assertion(!certainly(is_zero(sl)));
+
+      const typename K::FT t = ((v3v2.x()*v1v2.x()) + (v3v2.y()*v1v2.y()) + (v3v2.z()*v1v2.z())) / sl;
+      return t; // p1 + (p2-p1) * t
+    }
+  };
+
   template <typename K>
   class Construct_point_on_2
   {
diff --git a/Kernel_23/include/CGAL/Kernel/interface_macros.h b/Kernel_23/include/CGAL/Kernel/interface_macros.h
index 713e3773e56..a2314aef1b7 100644
--- a/Kernel_23/include/CGAL/Kernel/interface_macros.h
+++ b/Kernel_23/include/CGAL/Kernel/interface_macros.h
@@ -398,6 +398,12 @@ CGAL_Kernel_cons(Construct_plane_3,
                  construct_plane_3_object)
 CGAL_Kernel_cons(Construct_plane_line_intersection_point_3,
                  construct_plane_line_intersection_point_3_object)
+CGAL_Kernel_cons(Construct_planes_intersection_point_3,
+                 construct_planes_intersection_point_3_object)
+CGAL_Kernel_cons(Construct_coplanar_segments_intersection_point_3,
+                 construct_coplanar_segments_intersection_point_3_object)
+CGAL_Kernel_cons(Compute_alpha_for_coplanar_triangle_intersection_3,
+                 compute_alpha_for_coplanar_triangle_intersection_3_object)
 CGAL_Kernel_cons(Construct_point_on_2,
                  construct_point_on_2_object)
 CGAL_Kernel_cons(Construct_point_on_3,
diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Concepts/PMPAutorefinementVisitor.h b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Concepts/PMPAutorefinementVisitor.h
new file mode 100644
index 00000000000..ef2f124329b
--- /dev/null
+++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Concepts/PMPAutorefinementVisitor.h
@@ -0,0 +1,27 @@
+/// \ingroup PkgPolygonMeshProcessingConcepts
+/// \cgalConcept
+///
+/// The concept `PMPAutorefinementVisitor` defines the requirements for the visitor
+/// used in `CGAL::Polygon_mesh_processing::autorefine_triangle_soup()` to track
+/// the creation of new triangles.
+///
+/// \cgalRefines{CopyConstructible}
+/// \cgalHasModelsBegin
+/// \cgalHasModels{CGAL::Polygon_mesh_processing::Autorefinement::Default_visitor}
+/// \cgalHasModelsEnd
+
+class PMPAutorefinementVisitor{
+public:
+
+/// @name Functions called only if at least one intersection has been found
+/// @{
+  /// called when the final number of output triangles is known, `nbt` being the total number of triangles in the output.
+  void number_of_output_triangles(std::size_t nbt);
+  /// called for triangle with no intersection, `tgt_id` is the position in the triangle container after calling
+  /// `autorefine_triangle_soup()`, while `src_id` was its position before calling the function.
+  void verbatim_triangle_copy(std::size_t tgt_id, std::size_t src_id);
+  /// called for each subtriangle created from a triangle with intersection, `tgt_id` is the position in the triangle container after calling
+  /// `autorefine_triangle_soup()` of the subtriangle, while `src_id` was the position of the original support triangle before calling the function.
+  void new_subtriangle(std::size_t tgt_id, std::size_t src_id);
+/// @}
+};
diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt
index dc1da2b21fa..ad143bf835d 100644
--- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt
+++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt
@@ -113,6 +113,8 @@ The page \ref bgl_namedparameters "Named Parameters" describes their usage.
 - `CGAL::Polygon_mesh_processing::surface_intersection()`
 - `CGAL::Polygon_mesh_processing::clip()`
 - `CGAL::Polygon_mesh_processing::split()`
+- `CGAL::Polygon_mesh_processing::autorefine_triangle_soup()`
+- `CGAL::Polygon_mesh_processing::autorefine()`
 
 \cgalCRPSection{Meshing Functions}
 - `CGAL::Polygon_mesh_processing::remesh_planar_patches()`
diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt
index 26d5d9c6450..ac3e94db114 100644
--- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt
+++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt
@@ -229,7 +229,7 @@ it can try to equalize the angles between incident edges, or (and) move vertices
 Border vertices are considered constrained and do not move at any step of the procedure. No vertices are inserted at any time.
 Angle and area smoothing algorithms are based on Surazhsky and Gotsman \cgalCite{cgal:sg-hqct-04}.
 Since area smoothing considers only areas as a smoothing criterion, it may result in long and skinny
-triangles. To paliate this phenomenon, area smoothing is followed by an (optional) step of Delaunay-based edge flips.
+triangles. To palliate this phenomenon, area smoothing is followed by an (optional) step of Delaunay-based edge flips.
 In any case, area smoothing is guaranteed to improve the spatial distribution of the vertices
 over the area that is being smoothed. A simple example can be found in \ref Polygon_mesh_processing/mesh_smoothing_example.cpp.
 
@@ -449,7 +449,7 @@ depicted in red); (iii) `clip_volume=true`: clipping of the volume bounded by th
 \cgalFigureEnd
 
 \cgalFigureBegin{coref_clip_compact, clip_compact.png}
-Clipping a cube with a halfspace: compacity of the clipper (`clip_volume=false` in both cases).
+Clipping a cube with a halfspace: compactivity of the clipper (`clip_volume=false` in both cases).
 From left to right: (i) initial cube and the plane defining the clipping halfspace,
 note that a whole face of the cube (2 triangles) is exactly contained in the plane;
 (ii) `use_compact_clipper=true`: clipping of the surface mesh with a compact halfspace: coplanar faces are part of the output;
@@ -879,6 +879,19 @@ vertices at identical positions, can be used to repair this configuration.
 \section PMPGeometricRepair Geometric Repair
 ****************************************
 
+\subsection PMPAutoref Self-intersection Resolution (Autorefinement) in Triangle Soups
+Given a soup of triangles, a self-intersection is defined as the intersection of two triangles from the soup
+such that the intersection is not defined by the convex hull of one, two or three shared vertices.
+In other words, it is an intersection that happens in the interior of one of the two triangles, or in the interior
+of one their edges, except if identical points are associated to different vertices of the triangle soup which
+would then also includes overlaps of duplicated points.
+
+The function `CGAL::Polygon_mesh_processing::autorefine_triangle_soup()` provides a way to refine a triangle soup
+using the intersections of the triangles from the soup. In particular, if some points are duplicated they will be
+merged. Note that if a kernel with exact predicates but inexact constructions is used, some new self-intersections
+might be introduced due to rounding issues of points coordinates.
+To guarantee that the triangle soup is free from self-intersections, a kernel with exact constructions must be used.
+
 \subsection PMPRemoveCapsNeedles Removal of Almost Degenerate Triangle Faces
 Triangle faces of a mesh made up of almost collinear points are badly shaped elements that
 might not be desirable to have in a mesh. The function
diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt
index ba917e03836..3aca69fb004 100644
--- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt
+++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt
@@ -47,5 +47,6 @@
 \example Polygon_mesh_processing/remesh_planar_patches.cpp
 \example Polygon_mesh_processing/remesh_almost_planar_patches.cpp
 \example Polygon_mesh_processing/sample_example.cpp
+\example Polygon_mesh_processing/soup_autorefinement.cpp
 */
 */
diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt
index 90d9387d8ea..455e1ce2be7 100644
--- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt
+++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt
@@ -53,6 +53,8 @@ create_single_source_cgal_program("cc_compatible_orientations.cpp")
 create_single_source_cgal_program("hausdorff_distance_remeshing_example.cpp")
 create_single_source_cgal_program("hausdorff_bounded_error_distance_example.cpp")
 create_single_source_cgal_program("isotropic_remeshing_with_custom_sizing_example.cpp")
+create_single_source_cgal_program("triangle_mesh_autorefinement.cpp")
+create_single_source_cgal_program("soup_autorefinement.cpp")
 
 find_package(Eigen3 3.2.0 QUIET) #(requires 3.2.0 or greater)
 include(CGAL_Eigen3_support)
@@ -133,6 +135,7 @@ if(TARGET CGAL::TBB_support)
   target_link_libraries(self_intersections_example PUBLIC CGAL::TBB_support)
   target_link_libraries(hausdorff_distance_remeshing_example PUBLIC CGAL::TBB_support)
   target_link_libraries(hausdorff_bounded_error_distance_example PUBLIC CGAL::TBB_support)
+  target_link_libraries(soup_autorefinement PUBLIC CGAL::TBB_support)
 
   create_single_source_cgal_program("corefinement_parallel_union_meshes.cpp")
   target_link_libraries(corefinement_parallel_union_meshes PUBLIC CGAL::TBB_support)
diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/soup_autorefinement.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/soup_autorefinement.cpp
new file mode 100644
index 00000000000..a76e1254f6e
--- /dev/null
+++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/soup_autorefinement.cpp
@@ -0,0 +1,38 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Polygon_mesh_processing/autorefinement.h>
+#include <CGAL/Polygon_mesh_processing/repair_polygon_soup.h>
+#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
+#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
+#include <CGAL/IO/polygon_soup_io.h>
+#include <CGAL/Real_timer.h>
+
+#include <boost/container/small_vector.hpp>
+
+#include <iostream>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
+typedef Kernel::Point_3                                     Point;
+
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+int main(int argc, char** argv)
+{
+  const std::string filename = argc == 1 ? CGAL::data_file_path("meshes/elephant.off")
+                                         : std::string(argv[1]);
+
+  std::vector<Point> input_points;
+  std::vector<boost::container::small_vector<std::size_t, 3>> input_triangles;
+  CGAL::IO::read_polygon_soup(filename, input_points, input_triangles);
+  PMP::repair_polygon_soup(input_points, input_triangles);
+  PMP::triangulate_polygons(input_points, input_triangles);
+
+  CGAL::Real_timer t;
+  t.start();
+  PMP::autorefine_triangle_soup(input_points, input_triangles,
+                                CGAL::parameters::concurrency_tag(CGAL::Parallel_if_available_tag()));
+  t.stop();
+  std::cout << "#points = " << input_points.size() << " and #triangles = " << input_triangles.size() << " in " << t.time() << " sec." << std::endl;
+  CGAL::IO::write_polygon_soup("autorefined.off", input_points, input_triangles, CGAL::parameters::stream_precision(17));
+
+  return 0;
+}
diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/triangle_mesh_autorefinement.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/triangle_mesh_autorefinement.cpp
new file mode 100644
index 00000000000..65db52fef2c
--- /dev/null
+++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/triangle_mesh_autorefinement.cpp
@@ -0,0 +1,29 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Surface_mesh.h>
+
+#include <CGAL/Polygon_mesh_processing/autorefinement.h>
+#include <CGAL/boost/graph/IO/polygon_mesh_io.h>
+
+#include <iostream>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
+typedef Kernel::Point_3                                     Point;
+
+typedef CGAL::Surface_mesh<Point>                           Mesh;
+
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+int main(int argc, char** argv)
+{
+  Mesh mesh;
+
+  const std::string filename = argc == 1 ? CGAL::data_file_path("meshes/elephant.off")
+                                         : std::string(argv[1]);
+  CGAL::IO::read_polygon_mesh(filename, mesh);
+
+  PMP::autorefine(mesh);
+
+  CGAL::IO::write_polygon_mesh("autorefined.off", mesh, CGAL::parameters::stream_precision(17));
+
+  return 0;
+}
diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Non_manifold_feature_map.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Non_manifold_feature_map.h
index 8177a02d178..e8901ebd727 100644
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Non_manifold_feature_map.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Non_manifold_feature_map.h
@@ -28,8 +28,8 @@ struct Non_manifold_feature_map
   typedef typename Graph_traits::halfedge_descriptor halfedge_descriptor;
   typedef dynamic_edge_property_t<std::size_t> Edge_to_id_tag;
   typedef dynamic_vertex_property_t<std::size_t> Vertex_to_id_tag;
-  typedef typename boost::property_map<PolygonMesh, Edge_to_id_tag>::type Edge_to_nm_id;
-  typedef typename boost::property_map<PolygonMesh, Vertex_to_id_tag>::type Vertex_to_nm_id;
+  typedef typename boost::property_map<PolygonMesh, Edge_to_id_tag>::const_type Edge_to_nm_id;
+  typedef typename boost::property_map<PolygonMesh, Vertex_to_id_tag>::const_type Vertex_to_nm_id;
   Edge_to_nm_id e_nm_id;
   Vertex_to_nm_id v_nm_id;
   std::vector< std::vector<edge_descriptor> > non_manifold_edges;
@@ -39,7 +39,7 @@ struct Non_manifold_feature_map
   {}
 
   template <class Vpm>
-  Non_manifold_feature_map(PolygonMesh& pm, Vpm vpm)
+  Non_manifold_feature_map(const PolygonMesh& pm, Vpm vpm)
     : e_nm_id(get(Edge_to_id_tag(), pm))
     , v_nm_id(get(Vertex_to_id_tag(), pm))
   {
@@ -99,6 +99,14 @@ struct Non_manifold_feature_map
       }
     }
   }
+
+  void clear()
+  {
+    non_manifold_edges.clear();
+    non_manifold_vertices.clear();
+    e_nm_id = Edge_to_nm_id();
+    v_nm_id = Vertex_to_nm_id();
+  }
 };
 
 } } // end of CGAL::Polygon_mesh_processing
diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h
new file mode 100644
index 00000000000..8a61b92f717
--- /dev/null
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h
@@ -0,0 +1,1414 @@
+// Copyright (c) 2023 GeometryFactory (France).
+// 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)     : Sebastien Loriot
+//
+
+#ifndef CGAL_POLYGON_MESH_PROCESSING_AUTOREFINEMENT_H
+#define CGAL_POLYGON_MESH_PROCESSING_AUTOREFINEMENT_H
+
+#include <CGAL/license/Polygon_mesh_processing/geometric_repair.h>
+
+#include <CGAL/Polygon_mesh_processing/self_intersections.h>
+#include <CGAL/Polygon_mesh_processing/shape_predicates.h>
+#include <CGAL/Polygon_mesh_processing/repair_polygon_soup.h>
+#include <CGAL/Constrained_Delaunay_triangulation_2.h>
+#include <CGAL/Constrained_triangulation_plus_2.h>
+#include <CGAL/Projection_traits_3.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
+#ifdef CGAL_AUTOREF_USE_PROGRESS_DISPLAY
+#include <boost/timer/progress_display.hpp>
+#endif
+
+// output
+#include <CGAL/Polygon_mesh_processing/orient_polygon_soup_extension.h>
+#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
+#include <CGAL/Polygon_mesh_processing/polygon_mesh_to_polygon_soup.h>
+
+
+#ifdef CGAL_PMP_AUTOREFINE_USE_DEFAULT_VERBOSE
+#define CGAL_PMP_AUTOREFINE_VERBOSE(X) std::cout << X << "\n";
+#endif
+
+#ifndef CGAL_PMP_AUTOREFINE_VERBOSE
+#define CGAL_PMP_AUTOREFINE_VERBOSE(MSG)
+#endif
+
+#ifdef CGAL_LINKED_WITH_TBB
+#include <tbb/concurrent_vector.h>
+#if TBB_INTERFACE_VERSION < 12010 && !defined(TBB_PREVIEW_CONCURRENT_ORDERED_CONTAINERS)
+#define CGAL_HAS_DEFINED_TBB_PREVIEW_CONCURRENT_ORDERED_CONTAINERS
+#define TBB_PREVIEW_CONCURRENT_ORDERED_CONTAINERS 1
+#endif
+#include <tbb/concurrent_map.h>
+#include <tbb/parallel_for.h>
+#ifdef CGAL_AUTOREF_SET_POINT_IDS_USING_MUTEX
+#include <CGAL/mutex.h>
+#endif
+#endif
+
+#include <vector>
+
+//#define CGAL_AUTOREF_USE_DEBUG_PARALLEL_TIMERS
+//#define CGAL_AUTOREFINE_DEBUG_COUNTERS
+//#define CGAL_AUTOREF_DEBUG_DEPTH
+
+#ifdef CGAL_AUTOREF_USE_DEBUG_PARALLEL_TIMERS
+#include <CGAL/Real_timer.h>
+#endif
+
+#if defined(CGAL_AUTOREFINE_DEBUG_COUNTERS) || defined(CGAL_AUTOREF_USE_DEBUG_PARALLEL_TIMERS)
+#include <CGAL/Real_timer.h>
+#endif
+
+namespace CGAL {
+namespace Polygon_mesh_processing {
+
+namespace Autorefinement {
+
+/** \ingroup PMP_corefinement_grp
+ *  %Default visitor model of `PMPAutorefinementVisitor`.
+ *  All of its functions have an empty body. This class can be used as a
+ *  base class if only some of the functions of the concept require to be
+ *  overridden.
+ */
+struct Default_visitor
+{
+  inline void number_of_output_triangles(std::size_t /*nbt*/) {}
+  inline void verbatim_triangle_copy(std::size_t /*tgt_id*/, std::size_t /*src_id*/) {}
+  inline void new_subtriangle(std::size_t /*tgt_id*/, std::size_t /*src_id*/) {}
+};
+
+} // end of Autorefinement visitor
+
+
+#ifndef DOXYGEN_RUNNING
+namespace autorefine_impl {
+
+enum Segment_inter_type { NO_INTERSECTION=0,
+                          POINT_INTERSECTION,
+                          POINT_P,
+                          POINT_Q,
+                          POINT_R,
+                          POINT_S,
+                          COPLANAR_SEGMENT_PQ,
+                          COPLANAR_SEGMENT_RS,
+                          COPLANAR_SEGMENT_PS,
+                          COPLANAR_SEGMENT_QS,
+                          COPLANAR_SEGMENT_PR,
+                          COPLANAR_SEGMENT_QR,
+                        };
+
+// test intersection in the interior of segment pq and rs with pq and rs being coplanar segments
+// note that for coplanar cases, we might report identical endpoints
+template <class K>
+Segment_inter_type
+do_coplanar_segments_intersect(std::size_t pi, std::size_t qi,
+                               std::size_t ri, std::size_t si,
+                               const std::vector<typename K::Point_3>& points,
+                               const typename K::Vector_3& /* plane_normal */,
+                               const K& k = K())
+{
+  typename K::Collinear_are_ordered_along_line_3 cln_order = k.collinear_are_ordered_along_line_3_object();
+  typename K::Coplanar_orientation_3 cpl_orient=k.coplanar_orientation_3_object();
+
+  const typename K::Point_3& p=points[pi];
+  const typename K::Point_3& q=points[qi];
+  const typename K::Point_3& r=points[ri];
+  const typename K::Point_3& s=points[si];
+
+  // first handle case of shared endpoints
+  if (pi==ri)
+  {
+    if (si==qi || cpl_orient(p, q, s)!=COLLINEAR) return NO_INTERSECTION;
+    // can be s, q or nothing
+    if (cln_order(p,s,q))
+      return POINT_S;
+    if (cln_order(p,q,s))
+      return POINT_Q;
+    return NO_INTERSECTION;
+  }
+  else
+  {
+    if(pi==si)
+    {
+      if (qi==ri || cpl_orient(p, q, r)!=COLLINEAR) return NO_INTERSECTION;
+      // can be r, q or nothing
+      if (cln_order(p,r,q))
+        return POINT_R;
+      if (cln_order(p,q,r))
+        return POINT_Q;
+      return NO_INTERSECTION;
+    }
+    else
+    {
+      if (qi==ri)
+      {
+        if (pi==si || cpl_orient(p, q, s)!=COLLINEAR) return NO_INTERSECTION;
+        // can be p, s or nothing
+        if (cln_order(p,s,q))
+          return POINT_S;
+        if (cln_order(q,p,s))
+          return POINT_P;
+        return NO_INTERSECTION;
+      }
+      else
+      {
+        if (qi==si)
+        {
+          if (pi==ri || cpl_orient(p, q, r)!=COLLINEAR) return NO_INTERSECTION;
+          // can be p, r or nothing
+          if (cln_order(p,r,q))
+            return POINT_R;
+          if (cln_order(q,p,r))
+            return POINT_P;
+          return NO_INTERSECTION;
+        }
+      }
+    }
+  }
+
+  // supporting_line intersects: points are coplanar
+  ::CGAL::Orientation pqr = cpl_orient(p, q, r);
+  ::CGAL::Orientation pqs = cpl_orient(p, q, s);
+
+  if(pqr == COLLINEAR && pqs == COLLINEAR)
+  {
+    // segments are collinear
+    bool r_in_pq = cln_order(p, r, q),
+         s_in_pq = cln_order(p, s, q),
+         p_in_rs = cln_order(r, p, s);
+
+    if (r_in_pq)
+    {
+      // intersection could be rs, pr or qr
+      if (s_in_pq)
+        return COPLANAR_SEGMENT_RS;
+      if (p_in_rs)
+        return COPLANAR_SEGMENT_PR;
+      CGAL_assertion(cln_order(r, q, s));
+      return COPLANAR_SEGMENT_QR;
+    }
+    else
+    {
+      if (s_in_pq)
+      {
+        // intersection could be ps or qs
+        if (p_in_rs)
+          return COPLANAR_SEGMENT_PS;
+        CGAL_assertion(cln_order(r, q, s));
+        return COPLANAR_SEGMENT_QS;
+      }
+      else
+        if (p_in_rs)
+        {
+          CGAL_assertion(cln_order(r, q, s));
+          return COPLANAR_SEGMENT_PQ;
+        }
+    }
+    return NO_INTERSECTION;
+  }
+
+  if(pqr != pqs)
+  {
+    ::CGAL::Orientation rsp = cpl_orient(r, s, p);
+
+    if (rsp==COLLINEAR)
+    {
+      if (pqr==COLLINEAR || pqs==COLLINEAR)
+      {
+        throw std::runtime_error("no expected #1");
+      }
+      return POINT_P;
+    }
+    ::CGAL::Orientation rsq = cpl_orient(r, s, q);
+    if (rsq==COLLINEAR)
+    {
+      if (pqr==COLLINEAR || pqs==COLLINEAR)
+      {
+        throw std::runtime_error("no expected #2");
+      }
+      return POINT_Q;
+    }
+    if (rsp!=rsq)
+    {
+      if (pqr==COLLINEAR) return POINT_R;
+      if (pqs==COLLINEAR) return POINT_S;
+      return POINT_INTERSECTION;
+    }
+  }
+
+  return NO_INTERSECTION;
+}
+
+// imported from Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Triangle_3_intersection.h
+template<class K>
+void coplanar_intersections(const std::array<typename K::Point_3, 3>& t1,
+                            const std::array<typename K::Point_3, 3>& t2,
+                            std::vector<typename K::Point_3>& inter_pts)
+{
+  const typename K::Point_3& p1 = t1[0], q1 = t1[1], r1 = t1[2];
+  const typename K::Point_3& p2 = t2[0], q2 = t2[1], r2 = t2[2];
+
+  std::list<Intersections::internal::Point_on_triangle<K>> l_inter_pts;
+  l_inter_pts.push_back(Intersections::internal::Point_on_triangle<K>(-1,0));
+  l_inter_pts.push_back(Intersections::internal::Point_on_triangle<K>(-1,1));
+  l_inter_pts.push_back(Intersections::internal::Point_on_triangle<K>(-1,2));
+
+  //intersect t2 with the three half planes which intersection defines t1
+  K k;
+  Intersections::internal::intersection_coplanar_triangles_cutoff(p1,q1,r1,0,p2,q2,r2,k,l_inter_pts); //line p1q1
+  Intersections::internal::intersection_coplanar_triangles_cutoff(q1,r1,p1,1,p2,q2,r2,k,l_inter_pts); //line q1r1
+  Intersections::internal::intersection_coplanar_triangles_cutoff(r1,p1,q1,2,p2,q2,r2,k,l_inter_pts); //line r1p1
+
+  for (const Intersections::internal::Point_on_triangle<K>& pot : l_inter_pts)
+    inter_pts.push_back( pot.point(p1,q1,r1,p2,q2,r2,k) );
+}
+
+// imported from Polygon_mesh_processing/internal/Corefinement/intersect_triangle_segment_3.h
+template<class K>
+void
+find_intersection(const typename K::Point_3& p, const typename K::Point_3& q,  //segment
+                  const typename K::Point_3& a, const typename K::Point_3& b, const typename K::Point_3& c, //triangle
+                  std::vector<typename K::Point_3>& inter_pts,
+                  bool is_p_coplanar=false, bool is_q_coplanar=false) // note that in coref this was wrt a halfedge not p/q
+{
+  Orientation ab=orientation(p,q,a,b);
+  Orientation bc=orientation(p,q,b,c);
+  Orientation ca=orientation(p,q,c,a);
+
+  if ( ab==POSITIVE || bc==POSITIVE || ca==POSITIVE )
+    return;
+
+  int nb_coplanar=(ab==COPLANAR?1:0) + (bc==COPLANAR?1:0) + (ca==COPLANAR?1:0);
+
+  if (is_p_coplanar)
+  {
+    inter_pts.push_back(p);
+    return;
+  }
+  if (is_q_coplanar)
+  {
+    inter_pts.push_back(q);
+    return;
+  }
+
+  if (nb_coplanar!=2)
+  {
+    inter_pts.push_back(
+      typename K::Construct_plane_line_intersection_point_3()(a, b, c, p, q)
+    );
+  }
+  else
+  {
+    if (ab!=COPLANAR)
+    {
+      // intersection is c
+      inter_pts.push_back(c);
+      return;
+    }
+
+    if (bc!=COPLANAR)
+    {
+      // intersection is a
+      inter_pts.push_back(a);
+      return;
+    }
+    CGAL_assertion(ca!=COPLANAR);
+    // intersection is b
+    inter_pts.push_back(b);
+  }
+}
+
+template <class K>
+void test_edge(const typename K::Point_3& p, const typename K::Point_3& q,
+               const typename K::Point_3& a, const typename K::Point_3& b, const typename K::Point_3& c,
+               const Orientation abcp,
+               const Orientation abcq,
+               std::vector<typename K::Point_3>& inter_pts)
+{
+  switch ( abcp ) {
+  case POSITIVE:
+    switch ( abcq ) {
+    case POSITIVE:
+      // the segment lies in the positive open halfspaces defined by the
+      // triangle's supporting plane
+    break;
+    case NEGATIVE:
+      // p sees the triangle in counterclockwise order
+      find_intersection<K>(p,q,a,b,c,inter_pts);
+    break;
+    //case COPLANAR:
+    default:
+      // q belongs to the triangle's supporting plane
+      // p sees the triangle in counterclockwise order
+      find_intersection<K>(p,q,a,b,c,inter_pts,false,true);
+    }
+  break;
+  case NEGATIVE:
+    switch ( abcq ) {
+    case POSITIVE:
+      // q sees the triangle in counterclockwise order
+      find_intersection<K>(q,p,a,b,c,inter_pts);
+    break;
+    case NEGATIVE:
+      // the segment lies in the negative open halfspaces defined by the
+      // triangle's supporting plane
+    break;
+    // case COPLANAR:
+    default:
+      // q belongs to the triangle's supporting plane
+      // p sees the triangle in clockwise order
+      find_intersection<K>(q,p,a,b,c,inter_pts,true,false);
+    }
+  break;
+  default:
+  //case COPLANAR: // p belongs to the triangle's supporting plane
+    switch ( abcq ) {
+    case POSITIVE:
+      // q sees the triangle in counterclockwise order
+      find_intersection<K>(q,p,a,b,c,inter_pts,false, true);
+    break;
+    case NEGATIVE:
+      // q sees the triangle in clockwise order
+      find_intersection<K>(p,q,a,b,c,inter_pts,true);
+    break;
+    //case COPLANAR:
+    default:
+      // the segment is coplanar with the triangle's supporting plane
+      // we test whether the segment intersects the triangle in the common
+      // supporting plane
+      //if ( ::CGAL::Intersections::internal::do_intersect_coplanar(a,b,c,p,q,K()) )
+      //{
+        //handle coplanar intersection
+        // nothing done as coplanar case handle in collect_intersections
+        // and other intersection points will be collected with non-coplanar edges
+      //}
+    break;
+    }
+  }
+}
+
+template <class K>
+bool collect_intersections(const std::array<typename K::Point_3, 3>& t1,
+                           const std::array<typename K::Point_3, 3>& t2,
+                           std::vector<typename K::Point_3>& inter_pts)
+{
+  // test edges of t1 vs t2
+  std::array<Orientation,3> ori;
+  for (int i=0; i<3; ++i)
+    ori[i] = orientation(t2[0],t2[1],t2[2],t1[i]);
+
+  if (ori[0]== COPLANAR && ori[1]==COPLANAR && ori[2]==COPLANAR)
+  {
+    coplanar_intersections<K>(t1, t2, inter_pts);
+#ifdef CGAL_AUTOREF_DEBUG_DEPTH
+    for (auto p : inter_pts)
+      if (depth(p)>2) throw std::runtime_error("Depth is not 4: "+std::to_string(depth(p)));
+#endif
+
+    return true;
+  }
+
+  for (int i=0; i<3; ++i)
+  {
+    int j=(i+1)%3;
+    test_edge<K>(t1[i], t1[j], t2[0], t2[1], t2[2], ori[i], ori[j], inter_pts);
+  }
+
+  // test edges of t2 vs t1
+  for (int i=0; i<3; ++i)
+    ori[i] = orientation(t1[0],t1[1],t1[2],t2[i]);
+  for (int i=0; i<3; ++i)
+  {
+    int j=(i+1)%3;
+    test_edge<K>(t2[i], t2[j], t1[0], t1[1], t1[2], ori[i], ori[j], inter_pts);
+  }
+
+  // because we don't handle intersection type and can have edge-edge edge-vertex duplicates
+  std::sort(inter_pts.begin(), inter_pts.end());
+  auto last = std::unique(inter_pts.begin(), inter_pts.end());
+  inter_pts.erase(last, inter_pts.end());
+
+#ifdef CGAL_AUTOREF_DEBUG_DEPTH
+  for (auto p : inter_pts)
+    if (depth(p)>2) throw std::runtime_error("Depth is not 2: "+std::to_string(depth(p)));
+#endif
+
+  return false;
+}
+
+template <class EK,
+          class PointVector>
+void generate_subtriangles(std::size_t ti,
+                           std::vector<std::pair<std::size_t, std::size_t>>& segments,
+                           std::vector<typename EK::Point_3>& points,
+                           const std::vector<std::size_t>& in_triangle_ids,
+                           const std::set<std::pair<std::size_t, std::size_t> >& intersecting_triangles,
+                           const std::set<std::pair<std::size_t, std::size_t> >& coplanar_triangles,
+                           const std::vector<std::array<typename EK::Point_3,3>>& triangles,
+                           PointVector& new_triangles
+                           )
+{
+  typedef CGAL::Projection_traits_3<EK> P_traits;
+  typedef CGAL::No_constraint_intersection_tag Itag;
+
+  typedef CGAL::Constrained_Delaunay_triangulation_2<P_traits ,Default, Itag> CDT_2;
+  //typedef CGAL::Constrained_triangulation_plus_2<CDT_2> CDT;
+  typedef CDT_2 CDT;
+
+  const std::array<typename EK::Point_3,3>& t = triangles[ti];
+
+  // positive triangle normal
+  typename EK::Vector_3 n = normal(t[0], t[1], t[2]);
+  typename EK::Point_3 o(CGAL::ORIGIN);
+
+  bool orientation_flipped = false;
+  if ( typename EK::Less_xyz_3()(o+n,o) )
+  {
+    n=-n;
+    orientation_flipped = true;
+  }
+
+  P_traits cdt_traits(n);
+  CDT cdt(cdt_traits);
+  cdt.insert_outside_affine_hull(t[0]);
+  cdt.insert_outside_affine_hull(t[1]);
+  typename CDT::Vertex_handle v = cdt.tds().insert_dim_up(cdt.infinite_vertex(), orientation_flipped);
+  v->set_point(t[2]);
+
+#ifdef CGAL_AUTOREFINE_DEBUG_COUNTERS
+  struct Counter
+  {
+    int c1=0;
+    int c2=0;
+    int c3=0;
+    int c4=0;
+    int c5=0;
+    int total=0;
+    ~Counter()
+    {
+      std::cout << "intersection of 3 planes: " << c1 << "\n";
+      std::cout << "coplanar segment intersection : " << c2 << "\n";
+      std::cout << "coplanar segment overlap: " << c3 << "\n";
+      std::cout << "no intersection: " << c4 << "\n";
+      std::cout << "intersection filtered with bboxes: " << c5 << "\n";
+      std::cout << "# pairs of segments : " << total << "\n";
+      std::cout << "time computing segment intersections: " << timer1.time() << "\n";
+      std::cout << "time sorting intersection points: " << timer2.time() << "\n";
+      std::cout << "time for cdt of constraints: " << timer3.time() << "\n";
+      std::cout << "time coplanar segment intersections: " << timer4.time() << "\n";
+      std::cout << "time of do_coplanar_segments_intersect: " << timer5.time() << "\n";
+      std::cout << "time of triplane intersection: " << timer6.time() << "\n";
+    }
+    CGAL::Real_timer timer1, timer2, timer3, timer4, timer5, timer6;
+  };
+
+  static Counter counter;
+#define CGAL_AUTOREF_COUNTER_INSTRUCTION(X) X
+#else
+#define CGAL_AUTOREF_COUNTER_INSTRUCTION(X)
+#endif
+
+  // pre-compute segment intersections
+  if (!segments.empty())
+  {
+    std::size_t nbs = segments.size();
+
+    std::vector< std::vector<std::size_t> > points_on_segments(nbs);
+
+    CGAL_AUTOREF_COUNTER_INSTRUCTION(counter.timer1.start();)
+
+    std::map<typename EK::Point_3, std::size_t> point_id_map;
+
+    for (std::size_t pid=0; pid<points.size(); ++pid)
+    {
+      CGAL_assertion_code(bool insertion_ok =)
+      point_id_map.insert(std::make_pair(points[pid], pid))
+      CGAL_assertion_code(.second);
+      CGAL_assertion(insertion_ok);
+    }
+
+    auto get_point_id = [&](const typename EK::Point_3& pt)
+    {
+      auto insert_res = point_id_map.insert(std::make_pair(pt, points.size()));
+      if (insert_res.second)
+        points.push_back(pt);
+      return insert_res.first->second;
+    };
+
+    std::vector<Bbox_3> point_boxes(points.size());
+    for (std::size_t i = 0; i<points.size(); ++i)
+      point_boxes[i]=points[i].bbox();
+    std::vector<Bbox_3> segment_boxes(nbs);
+    for (std::size_t i = 0; i<nbs; ++i)
+      segment_boxes[i]=point_boxes[segments[i].first]+point_boxes[segments[i].second];
+
+    for (std::size_t i = 0; i<nbs-1; ++i)
+    {
+      for (std::size_t j = i+1; j<nbs; ++j)
+      {
+        if (intersecting_triangles.count(CGAL::make_sorted_pair(in_triangle_ids[i], in_triangle_ids[j]))!=0)
+        {
+          if ( !do_overlap(segment_boxes[i], segment_boxes[j]) )
+          {
+            CGAL_AUTOREF_COUNTER_INSTRUCTION(++counter.c5;)
+            CGAL_AUTOREF_COUNTER_INSTRUCTION(++counter.total;)
+            continue;
+          }
+
+          CGAL_AUTOREF_COUNTER_INSTRUCTION(counter.timer5.start();)
+          Segment_inter_type seg_inter_type =
+            do_coplanar_segments_intersect<EK>(segments[i].first, segments[i].second,
+                                               segments[j].first, segments[j].second,
+                                               points, n);
+          CGAL_AUTOREF_COUNTER_INSTRUCTION(counter.timer5.stop();)
+
+          switch(seg_inter_type)
+          {
+            case POINT_P:
+            {
+              points_on_segments[j].push_back(segments[i].first);
+              break;
+            }
+            case POINT_Q:
+            {
+              points_on_segments[j].push_back(segments[i].second);
+              break;
+            }
+            case POINT_R:
+            {
+              points_on_segments[i].push_back(segments[j].first);
+              break;
+            }
+            case POINT_S:
+            {
+              points_on_segments[i].push_back(segments[j].second);
+              break;
+            }
+            case POINT_INTERSECTION:
+            {
+              if (   coplanar_triangles.count(CGAL::make_sorted_pair(in_triangle_ids[i], in_triangle_ids[j])) == 0
+                  && coplanar_triangles.count(CGAL::make_sorted_pair(ti, in_triangle_ids[j])) == 0
+                  && coplanar_triangles.count(CGAL::make_sorted_pair(in_triangle_ids[i], ti)) == 0)
+              {
+                CGAL_AUTOREF_COUNTER_INSTRUCTION(counter.timer6.start();)
+                typename EK::Point_3 pt = typename EK::Construct_planes_intersection_point_3()(
+                  triangles[in_triangle_ids[i]][0], triangles[in_triangle_ids[i]][1],triangles[in_triangle_ids[i]][2],
+                  triangles[in_triangle_ids[j]][0], triangles[in_triangle_ids[j]][1],triangles[in_triangle_ids[j]][2],
+                  triangles[ti][0], triangles[ti][1],triangles[ti][2]);
+                CGAL_AUTOREF_COUNTER_INSTRUCTION(counter.timer6.stop();)
+
+                CGAL_AUTOREF_COUNTER_INSTRUCTION(++counter.c1;)
+                std::size_t pid = get_point_id(pt);
+                points_on_segments[i].push_back(pid);
+                points_on_segments[j].push_back(pid);
+              }
+              else
+              {
+                CGAL_AUTOREF_COUNTER_INSTRUCTION(++counter.c2;)
+                CGAL_AUTOREF_COUNTER_INSTRUCTION(counter.timer4.start();)
+                typename EK::Point_3 pt = typename EK::Construct_coplanar_segments_intersection_point_3()(
+                  points[segments[i].first], points[segments[i].second],
+                  points[segments[j].first], points[segments[j].second]);
+
+                std::size_t pid = get_point_id(pt);
+                points_on_segments[i].push_back(pid);
+                points_on_segments[j].push_back(pid);
+                CGAL_AUTOREF_COUNTER_INSTRUCTION(counter.timer4.stop();)
+              }
+              break;
+            }
+            case COPLANAR_SEGMENT_PQ:
+            {
+              CGAL_AUTOREF_COUNTER_INSTRUCTION(++counter.c3;)
+              points_on_segments[j].push_back(segments[i].first);
+              points_on_segments[j].push_back(segments[i].second);
+              break;
+            }
+            case COPLANAR_SEGMENT_RS:
+            {
+              CGAL_AUTOREF_COUNTER_INSTRUCTION(++counter.c3;)
+              points_on_segments[i].push_back(segments[j].first);
+              points_on_segments[i].push_back(segments[j].second);
+              break;
+            }
+            case COPLANAR_SEGMENT_PR:
+            {
+              CGAL_AUTOREF_COUNTER_INSTRUCTION(++counter.c3;)
+              points_on_segments[i].push_back(segments[j].first);
+              points_on_segments[j].push_back(segments[i].first);
+              break;
+            }
+            case COPLANAR_SEGMENT_QS:
+            {
+              CGAL_AUTOREF_COUNTER_INSTRUCTION(++counter.c3;)
+              points_on_segments[i].push_back(segments[j].second);
+              points_on_segments[j].push_back(segments[i].second);
+              break;
+            }
+            case COPLANAR_SEGMENT_PS:
+            {
+              CGAL_AUTOREF_COUNTER_INSTRUCTION(++counter.c3;)
+              points_on_segments[i].push_back(segments[j].second);
+              points_on_segments[j].push_back(segments[i].first);
+              break;
+            }
+            case COPLANAR_SEGMENT_QR:
+            {
+              CGAL_AUTOREF_COUNTER_INSTRUCTION(++counter.c3;)
+              points_on_segments[i].push_back(segments[j].first);
+              points_on_segments[j].push_back(segments[i].second);
+              break;
+            }
+            case NO_INTERSECTION:
+            {
+              CGAL_AUTOREF_COUNTER_INSTRUCTION(++counter.c4;)
+            }
+          }
+        }
+        CGAL_AUTOREF_COUNTER_INSTRUCTION(++counter.total;)
+      }
+    }
+    CGAL_AUTOREF_COUNTER_INSTRUCTION(counter.timer1.stop();)
+    CGAL_AUTOREF_COUNTER_INSTRUCTION(counter.timer2.start();)
+    std::size_t nb_new_segments=0;
+    for (std::size_t i = 0; i<nbs; ++i)
+    {
+      if(!points_on_segments[i].empty())
+      {
+        int coord = 0;
+        std::size_t src_id = segments[i].first, tgt_id = segments[i].second;
+        typename EK::Point_3 src = points[src_id], tgt=points[tgt_id];
+        if (src.x()==tgt.x())
+        {
+          coord=1;
+          if (src.y()==tgt.y())
+            coord=2;
+        }
+        if (src[coord]>tgt[coord])
+        {
+          std::swap(src_id, tgt_id);
+          std::swap(src, tgt);
+        }
+
+        points_on_segments[i].push_back(src_id);
+        std::swap(points_on_segments[i].front(), points_on_segments[i].back());
+        std::sort(std::next(points_on_segments[i].begin()), points_on_segments[i].end(),
+                  [&](std::size_t id1, std::size_t id2)
+                  {
+                    if (id1==id2) return false;
+                    return points[id1][coord]<points[id2][coord];
+                  });
+        points_on_segments[i].push_back(tgt_id);
+        auto last = std::unique(points_on_segments[i].begin(), points_on_segments[i].end());
+        points_on_segments[i].erase(last, points_on_segments[i].end());
+        nb_new_segments+=points_on_segments[i].size()-2;
+      }
+    }
+    CGAL_AUTOREF_COUNTER_INSTRUCTION(counter.timer2.stop();)
+
+    // now fill segments with new segments
+    segments.reserve(segments.size()+nb_new_segments);
+    for (std::size_t i = 0; i<nbs; ++i)
+    {
+      if(!points_on_segments[i].empty())
+      {
+        segments[i]=std::make_pair(points_on_segments[i][0], points_on_segments[i][1]);
+        for(std::size_t pos=1; pos<points_on_segments[i].size()-1; ++pos)
+          segments.emplace_back(points_on_segments[i][pos], points_on_segments[i][pos+1]);
+      }
+    }
+
+    CGAL_assertion(points.size()==std::set<typename EK::Point_3>(points.begin(), points.end()).size());
+    CGAL_assertion(points.size()==point_id_map.size());
+  }
+
+  for (std::pair<std::size_t, size_t>& s : segments)
+    if (s.second < s.first)
+      std::swap(s.first,s.second);
+  std::sort(segments.begin(), segments.end());
+  auto last = std::unique(segments.begin(), segments.end());
+  segments.erase(last, segments.end());
+
+  CGAL_AUTOREF_COUNTER_INSTRUCTION(counter.timer3.start();)
+  if (segments.empty())
+    cdt.insert(points.begin(), points.end());
+  else
+    cdt.insert_constraints(points.begin(), points.end(), segments.begin(), segments.end());
+  CGAL_AUTOREF_COUNTER_INSTRUCTION(counter.timer3.stop();)
+
+  for (typename CDT::Face_handle fh : cdt.finite_face_handles())
+  {
+    if (orientation_flipped)
+      new_triangles.push_back( { CGAL::make_array(fh->vertex(0)->point(),
+                                                  fh->vertex(cdt.cw(0))->point(),
+                                                  fh->vertex(cdt.ccw(0))->point()), ti } );
+    else
+      new_triangles.push_back( { CGAL::make_array(fh->vertex(0)->point(),
+                                                  fh->vertex(cdt.ccw(0))->point(),
+                                                  fh->vertex(cdt.cw(0))->point()), ti } );
+  }
+}
+
+} // end of autorefine_impl
+#endif
+
+/**
+* \ingroup PMP_corefinement_grp
+*
+* refines a soup of triangles so that no pair of triangles intersects.
+* Output triangles may share a common edge or a common vertex (but with the same indexed position in `points`).
+* Note that points in `soup_points` can only be added (intersection points) at the end of the container, with the initial order preserved.
+* Note that if `soup_points` contains two or more identical points then only the first copy (following the order in the `soup_points`)
+* will be used in `soup_triangles`.
+* `soup_triangles` will be updated to contain both the input triangles and the new subdivided triangles. Degenerate triangles will be removed.
+* Also triangles in `soup_triangles` will be triangles without intersection first, followed by triangles coming from a subdivision induced
+* by an intersection. The named parameter `visitor()` can be used to track
+*
+* @tparam PointRange a model of the concept `RandomAccessContainer`
+* whose value type is the point type
+* @tparam TriangleRange a model of the concepts `RandomAccessContainer`, `BackInsertionSequence` and `Swappable`, whose
+* value type is a model of the concept `RandomAccessContainer` whose value type is convertible to `std::size_t` and that
+* is constructible from an `std::initializer_list<std::size_t>` of size 3.
+* @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
+*
+* @param soup_points points of the soup of polygons
+* @param soup_triangles each element in the range describes a triangle using the indexed position of the points in `soup_points`
+* @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+*
+* \cgalNamedParamsBegin
+*   \cgalParamNBegin{concurrency_tag}
+*     \cgalParamDescription{a tag indicating if the task should be done using one or several threads.}
+*     \cgalParamType{Either `CGAL::Sequential_tag`, or `CGAL::Parallel_tag`, or `CGAL::Parallel_if_available_tag`}
+*     \cgalParamDefault{`CGAL::Sequential_tag`}
+*   \cgalParamNEnd
+*   \cgalParamNBegin{point_map}
+*     \cgalParamDescription{a property map associating points to the elements of the range `soup_points`}
+*     \cgalParamType{a model of `ReadWritePropertyMap` whose value type is a point type}
+*     \cgalParamDefault{`CGAL::Identity_property_map`}
+*   \cgalParamNEnd
+*   \cgalParamNBegin{geom_traits}
+*     \cgalParamDescription{an instance of a geometric traits class}
+*     \cgalParamType{a class model of `Kernel`}
+*     \cgalParamDefault{a \cgal Kernel deduced from the point type, using `CGAL::Kernel_traits`}
+*     \cgalParamExtra{The geometric traits class must be compatible with the point type.}
+*   \cgalParamNEnd
+*   \cgalParamNBegin{visitor}
+*     \cgalParamDescription{a visitor used to track the creation of new faces}
+*     \cgalParamType{a class model of `PMPAutorefinementVisitor`}
+*     \cgalParamDefault{`Autorefinement::Default_visitor`}
+*     \cgalParamExtra{The visitor will be copied.}
+*   \cgalParamNEnd
+* \cgalNamedParamsEnd
+*
+*/
+template <class PointRange, class TriangleRange, class NamedParameters = parameters::Default_named_parameters>
+void autorefine_triangle_soup(PointRange& soup_points,
+                              TriangleRange& soup_triangles,
+                              const NamedParameters& np = parameters::default_values())
+{
+  using parameters::choose_parameter;
+  using parameters::get_parameter;
+
+  typedef typename GetPolygonSoupGeomTraits<PointRange, NamedParameters>::type GT;
+  typedef typename GetPointMap<PointRange, NamedParameters>::const_type    Point_map;
+  Point_map pm = choose_parameter<Point_map>(get_parameter(np, internal_np::point_map));
+
+  typedef typename internal_np::Lookup_named_param_def <
+    internal_np::concurrency_tag_t,
+    NamedParameters,
+    Sequential_tag
+  > ::type Concurrency_tag;
+
+  // visitor
+  typedef typename internal_np::Lookup_named_param_def <
+    internal_np::visitor_t,
+    NamedParameters,
+    Autorefinement::Default_visitor//default
+  > ::type Visitor;
+  Visitor visitor(choose_parameter<Visitor>(get_parameter(np, internal_np::visitor)));
+
+
+  constexpr bool parallel_execution = std::is_same_v<Parallel_tag, Concurrency_tag>;
+
+#ifndef CGAL_LINKED_WITH_TBB
+  static_assert (!parallel_execution,
+                 "Parallel_tag is enabled but TBB is unavailable.");
+#endif
+
+  typedef std::size_t Input_TID;
+  typedef std::pair<Input_TID, Input_TID> Pair_of_triangle_ids;
+
+  // no need for a concurrent vector as the called function itself
+  // takes care of sequentially writing into the output iterator
+  std::vector<Pair_of_triangle_ids> si_pairs;
+
+  // collect intersecting pairs of triangles
+  CGAL_PMP_AUTOREFINE_VERBOSE("collect intersecting pairs");
+  triangle_soup_self_intersections<Concurrency_tag>(soup_points, soup_triangles, std::back_inserter(si_pairs), np);
+
+  if (si_pairs.empty())
+  {
+    if constexpr (!std::is_same_v<Autorefinement::Default_visitor, Visitor>)
+    {
+      visitor.number_of_output_triangles(soup_triangles.size());
+      for(std::size_t i=0; i<soup_triangles.size(); ++i)
+        visitor.verbatim_triangle_copy(i, i);
+    }
+    return;
+  }
+
+  // mark degenerate faces so that we can ignore them
+  std::vector<bool> is_degen(soup_triangles.size(), false);
+
+  for (const Pair_of_triangle_ids& p : si_pairs)
+    if (p.first==p.second) // bbox inter reports (f,f) for degenerate faces
+      is_degen[p.first] = true;
+
+  // assign an id per triangle involved in an intersection
+  // + the faces involved in the intersection
+  std::vector<int> tri_inter_ids(soup_triangles.size(), -1);
+  std::vector<Input_TID> intersected_faces;
+  int tiid=-1;
+  for (const Pair_of_triangle_ids& p : si_pairs)
+  {
+    if (tri_inter_ids[p.first]==-1 && !is_degen[p.first])
+    {
+      tri_inter_ids[p.first]=++tiid;
+      intersected_faces.push_back(p.first);
+    }
+    if (tri_inter_ids[p.second]==-1 && !is_degen[p.second])
+    {
+      tri_inter_ids[p.second]=++tiid;
+      intersected_faces.push_back(p.second);
+    }
+  }
+
+  // init the vector of triangles used for the autorefinement of triangles
+  typedef CGAL::Exact_predicates_exact_constructions_kernel EK;
+  std::vector< std::array<EK::Point_3, 3> > triangles(tiid+1);
+  Cartesian_converter<GT, EK> to_exact;
+
+  for(Input_TID f : intersected_faces)
+  {
+    triangles[tri_inter_ids[f]]= CGAL::make_array(
+      to_exact( get(pm, soup_points[soup_triangles[f][0]]) ),
+      to_exact( get(pm, soup_points[soup_triangles[f][1]]) ),
+      to_exact( get(pm, soup_points[soup_triangles[f][2]]) ) );
+  }
+
+  std::vector< std::vector<std::array<EK::Point_3,2> > > all_segments(triangles.size());
+  std::vector< std::vector<EK::Point_3> > all_points(triangles.size());
+  std::vector< std::vector<std::size_t> > all_in_triangle_ids(triangles.size());
+
+  CGAL_PMP_AUTOREFINE_VERBOSE("compute intersections");
+#ifdef CGAL_AUTOREF_USE_DEBUG_PARALLEL_TIMERS
+  Real_timer t;
+  t.start();
+#endif
+  std::set<std::pair<std::size_t, std::size_t> > intersecting_triangles;
+  std::set<std::pair<std::size_t, std::size_t> > coplanar_triangles;
+  //TODO: PARALLEL_FOR #2
+  for (const Pair_of_triangle_ids& p : si_pairs)
+  {
+    int i1 = tri_inter_ids[p.first],
+        i2 = tri_inter_ids[p.second];
+
+    if (i1==-1 || i2==-1) continue; //skip degenerate faces
+
+    const std::array<EK::Point_3, 3>& t1 = triangles[i1];
+    const std::array<EK::Point_3, 3>& t2 = triangles[i2];
+
+    std::vector<EK::Point_3> inter_pts;
+    bool triangles_are_coplanar = autorefine_impl::collect_intersections<EK>(t1, t2, inter_pts);
+
+    CGAL_assertion(
+      CGAL::do_intersect(EK::Triangle_3(t1[0], t1[1], t1[2]), EK::Triangle_3(t2[0], t2[1], t2[2]))
+      != inter_pts.empty());
+
+    if (!inter_pts.empty())
+    {
+      std::size_t nbi = inter_pts.size();
+      switch(nbi)
+      {
+        case 1:
+          all_points[i1].push_back(inter_pts[0]);
+          all_points[i2].push_back(inter_pts[0]);
+        break;
+        case 2:
+          all_segments[i1].push_back(CGAL::make_array(inter_pts[0], inter_pts[1]));
+          all_segments[i2].push_back(CGAL::make_array(inter_pts[0], inter_pts[1]));
+          all_in_triangle_ids[i1].push_back(i2);
+          all_in_triangle_ids[i2].push_back(i1);
+        break;
+        default:
+          for (std::size_t i=0;i<nbi; ++i)
+          {
+            all_segments[i1].push_back(CGAL::make_array(inter_pts[i], inter_pts[(i+1)%nbi]));
+            all_segments[i2].push_back(all_segments[i1].back());
+            all_in_triangle_ids[i1].push_back(i2);
+            all_in_triangle_ids[i2].push_back(i1);
+          }
+      }
+      intersecting_triangles.insert(CGAL::make_sorted_pair(i1, i2));
+      if (triangles_are_coplanar)
+        coplanar_triangles.insert(CGAL::make_sorted_pair(i1, i2));
+    }
+  }
+#ifdef CGAL_AUTOREF_USE_DEBUG_PARALLEL_TIMERS
+  t.stop();
+  std::cout << t.time() << " sec. for #2" << std::endl;
+  t.reset();
+#endif
+
+#ifdef CGAL_AUTOREF_USE_DEBUG_PARALLEL_TIMERS
+  t.start();
+#endif
+
+  // deduplicate inserted segments
+  std::vector<std::vector<std::pair<std::size_t, std::size_t>>> all_segments_ids(all_segments.size());
+
+  auto deduplicate_inserted_segments = [&](std::size_t ti)
+  {
+    if (!all_segments[ti].empty())
+    {
+      std::map<EK::Point_3, std::size_t> point_id_map;
+
+
+      auto get_point_id = [&](const EK::Point_3& pt)
+      {
+        auto insert_res = point_id_map.insert(std::make_pair(pt, all_points[ti].size()));
+        if (insert_res.second)
+          all_points[ti].push_back(pt);
+        return insert_res.first->second;
+      };
+
+
+      if (!all_points[ti].empty())
+      {
+        using EPoint_3 = EK::Point_3; // workaround for MSVC 2022 bug
+        std::vector<EPoint_3> tmp;
+        tmp.swap(all_points[ti]);
+        for (const EPoint_3& pt : tmp)
+          get_point_id(pt);
+      }
+
+      std::size_t nbs = all_segments[ti].size();
+      std::vector<std::array<EK::Point_3,2>> filtered_segments;
+      std::vector<std::size_t> filtered_in_triangle_ids;
+      filtered_segments.reserve(nbs);
+      std::set<std::pair<std::size_t, std::size_t>> segset;
+      for (std::size_t si=0; si<nbs; ++si)
+      {
+        const EK::Point_3& src = all_segments[ti][si][0];
+        const EK::Point_3& tgt = all_segments[ti][si][1];
+        std::size_t src_id = get_point_id(src), tgt_id=get_point_id(tgt);
+        if (segset.insert(
+              CGAL::make_sorted_pair(src_id, tgt_id)).second)
+        {
+          all_segments_ids[ti].emplace_back(src_id, tgt_id);
+          filtered_in_triangle_ids.push_back(all_in_triangle_ids[ti][si]);
+        }
+      }
+      if (all_segments_ids[ti].size()!=nbs)
+        filtered_in_triangle_ids.swap(all_in_triangle_ids[ti]);
+    }
+  };
+
+#ifdef CGAL_LINKED_WITH_TBB
+  if (parallel_execution)
+  {
+    tbb::parallel_for(tbb::blocked_range<size_t>(0, triangles.size()),
+                      [&](const tbb::blocked_range<size_t>& r) {
+                        for (size_t ti = r.begin(); ti != r.end(); ++ti)
+                          deduplicate_inserted_segments(ti);
+                      }
+                      );
+  }
+  else
+#endif
+  for (std::size_t ti = 0; ti < triangles.size(); ++ti) {
+    deduplicate_inserted_segments(ti);
+  }
+
+#ifdef CGAL_AUTOREF_USE_DEBUG_PARALLEL_TIMERS
+  t.stop();
+  std::cout << t.time() << " sec. for #3" << std::endl;
+  t.reset();
+#endif
+
+  CGAL_PMP_AUTOREFINE_VERBOSE("triangulate faces");
+  // now refine triangles
+#ifdef CGAL_LINKED_WITH_TBB
+  std::conditional_t<parallel_execution,
+                     tbb::concurrent_vector<std::pair<std::array<EK::Point_3,3>, std::size_t>>,
+                     std::vector<std::pair<std::array<EK::Point_3,3>, std::size_t>>> new_triangles;
+#else
+  std::vector<std::pair<std::array<EK::Point_3,3>, std::size_t>> new_triangles;
+#endif
+
+#ifdef CGAL_AUTOREF_USE_PROGRESS_DISPLAY
+  boost::timer::progress_display pd(triangles.size());
+#endif
+
+  auto refine_triangles = [&](std::size_t ti)
+  {
+    if (all_segments[ti].empty() && all_points[ti].empty())
+      new_triangles.push_back({triangles[ti], ti});
+    else
+    {
+      autorefine_impl::generate_subtriangles<EK>(ti, all_segments_ids[ti], all_points[ti], all_in_triangle_ids[ti], intersecting_triangles, coplanar_triangles, triangles, new_triangles);
+    }
+
+#ifdef CGAL_AUTOREF_USE_PROGRESS_DISPLAY
+    ++pd;
+#endif
+  };
+
+#ifdef CGAL_AUTOREF_USE_DEBUG_PARALLEL_TIMERS
+  t.start();
+#endif
+#ifdef CGAL_LINKED_WITH_TBB
+  if (parallel_execution)
+  {
+    tbb::parallel_for(tbb::blocked_range<size_t>(0, triangles.size()),
+                      [&](const tbb::blocked_range<size_t>& r) {
+                        for (size_t ti = r.begin(); ti != r.end(); ++ti)
+                          refine_triangles(ti);
+                      }
+                      );
+  }
+  else
+#endif
+  for (std::size_t ti = 0; ti < triangles.size(); ++ti) {
+    refine_triangles(ti);
+  }
+
+#ifdef CGAL_AUTOREF_USE_DEBUG_PARALLEL_TIMERS
+  t.stop();
+  std::cout << t.time() << " sec. for #1" << std::endl;
+  t.reset();
+#endif
+
+  // brute force output: create a soup, orient and to-mesh
+  CGAL_PMP_AUTOREFINE_VERBOSE("create output soup");
+
+  Cartesian_converter<EK, GT> to_input;
+
+#ifdef CGAL_LINKED_WITH_TBB
+  typedef std::conditional_t<parallel_execution,
+                             tbb::concurrent_map<EK::Point_3, std::size_t>,
+                             std::map<EK::Point_3, std::size_t>> Point_id_map;
+#else
+  typedef std::map<EK::Point_3, std::size_t> Point_id_map;
+#endif
+  Point_id_map point_id_map;
+
+#if ! defined(CGAL_NDEBUG) || defined(CGAL_DEBUG_PMP_AUTOREFINE)
+  std::vector<EK::Point_3> exact_soup_points;
+#endif
+
+  // TODO: parallel_for?
+  // for input points, we on purpose keep duplicated points and isolated points
+  for (std::size_t pid = 0; pid<soup_points.size(); ++pid)
+  {
+#if ! defined(CGAL_NDEBUG) || defined(CGAL_DEBUG_PMP_AUTOREFINE)
+    auto insert_res =
+#endif
+    point_id_map.insert(
+      std::make_pair(to_exact(get(pm,soup_points[pid])), pid));
+#if ! defined(CGAL_NDEBUG) || defined(CGAL_DEBUG_PMP_AUTOREFINE)
+      exact_soup_points.push_back(insert_res.first->first);
+#endif
+  }
+
+  TriangleRange soup_triangles_out;
+  soup_triangles_out.reserve(soup_triangles.size());
+
+  if constexpr (!std::is_same_v<Autorefinement::Default_visitor, Visitor>)
+  {
+    std::size_t nbt=0;
+    for (Input_TID f=0; f<soup_triangles.size(); ++f)
+    {
+      int tiid = tri_inter_ids[f];
+      if (tiid == -1) ++nbt;
+    }
+    visitor.number_of_output_triangles(nbt+new_triangles.size());
+  }
+
+  // raw copy of input triangles with no intersection
+  std::vector<std::size_t> tri_inter_ids_inverse(triangles.size());
+  for (Input_TID f=0; f<soup_triangles.size(); ++f)
+  {
+    if (is_degen[f]) continue; //skip degenerate faces
+
+    int tiid = tri_inter_ids[f];
+    if (tiid == -1)
+    {
+      visitor.verbatim_triangle_copy(soup_triangles_out.size(), f);
+      soup_triangles_out.push_back(
+        {soup_triangles[f][0], soup_triangles[f][1], soup_triangles[f][2]}
+      );
+    }
+    else
+    {
+      tri_inter_ids_inverse[tiid]=f;
+    }
+  }
+
+  // import refined triangles
+
+  // Lambda to retrieve the id of intersection points
+  auto get_point_id = [&](const EK::Point_3& pt)
+  {
+    auto insert_res = point_id_map.insert(std::make_pair(pt, soup_points.size()));
+    if (insert_res.second)
+    {
+      soup_points.push_back(to_input(pt));
+#if ! defined(CGAL_NDEBUG) || defined(CGAL_DEBUG_PMP_AUTOREFINE)
+      exact_soup_points.push_back(pt);
+#endif
+    }
+    return insert_res.first->second;
+  };
+
+#ifdef CGAL_AUTOREF_USE_DEBUG_PARALLEL_TIMERS
+  t.start();
+#endif
+
+  std::size_t offset = soup_triangles_out.size();
+#ifdef CGAL_AUTOREF_USE_DEBUG_PARALLEL_TIMERS
+  std::string mode = "parallel";
+#endif
+
+// It might be possible to optimise the hardcoded value below
+// but the less triangles the faster will anyway be the operation.
+// So it's probably not critical.
+#ifdef CGAL_LINKED_WITH_TBB
+  if(parallel_execution && new_triangles.size() > 50)
+  {
+#ifdef CGAL_AUTOREF_SET_POINT_IDS_USING_MUTEX
+    //option 1 (using a mutex)
+    CGAL_MUTEX point_container_mutex;
+    /// Lambda concurrent_get_point_id()
+    auto concurrent_get_point_id = [&](const EK::Point_3& pt)
+    {
+      auto insert_res = point_id_map.insert(std::make_pair(pt, -1));
+
+      if (insert_res.second)
+      {
+        CGAL_SCOPED_LOCK(point_container_mutex);
+        insert_res.first->second=soup_points.size();
+        soup_points.push_back(to_input(pt));
+#if ! defined(CGAL_NDEBUG) || defined(CGAL_DEBUG_PMP_AUTOREFINE)
+        exact_soup_points.push_back(pt);
+#endif
+      }
+      return insert_res.first;
+    };
+
+    soup_triangles_out.resize(offset + new_triangles.size());
+    //use map iterator triple for triangles to create them concurrently and safely
+    std::vector<std::array<typename Point_id_map::iterator, 3>> triangle_buffer(new_triangles.size());
+    tbb::parallel_for(tbb::blocked_range<size_t>(0, new_triangles.size()),
+      [&](const tbb::blocked_range<size_t>& r) {
+        for (size_t ti = r.begin(); ti != r.end(); ++ti) {
+          const std::array<EK::Point_3, 3>& t = new_triangles[ti].first;
+          visitor.new_subtriangle(offset+ti, tri_inter_ids_inverse[new_triangles[ti].second]);
+          triangle_buffer[ti] = CGAL::make_array(concurrent_get_point_id(t[0]), concurrent_get_point_id(t[1]), concurrent_get_point_id(t[2]));
+        }
+      }
+    );
+    tbb::parallel_for(tbb::blocked_range<size_t>(0, new_triangles.size()),
+      [&](const tbb::blocked_range<size_t>& r) {
+          for (size_t ti = r.begin(); ti != r.end(); ++ti)
+          {
+            soup_triangles_out[offset + ti] =
+              { triangle_buffer[ti][0]->second,
+                triangle_buffer[ti][1]->second,
+                triangle_buffer[ti][2]->second };
+          }
+        }
+    );
+#else
+    //option 2 (without mutex)
+    /// Lambda concurrent_get_point_id()
+    tbb::concurrent_vector<typename Point_id_map::iterator> iterators;
+    auto concurrent_get_point_id = [&](const EK::Point_3& pt)
+    {
+      auto insert_res = point_id_map.insert(std::make_pair(pt, -1));
+      if (insert_res.second)
+        iterators.push_back(insert_res.first);
+      return insert_res.first;
+    };
+
+    //use map iterator triple for triangles to create them concurrently and safely
+    soup_triangles_out.resize(offset + new_triangles.size());
+    std::vector<std::array<typename Point_id_map::iterator, 3>> triangle_buffer(new_triangles.size());
+    tbb::parallel_for(tbb::blocked_range<size_t>(0, new_triangles.size()),
+      [&](const tbb::blocked_range<size_t>& r) {
+        for (size_t ti = r.begin(); ti != r.end(); ++ti)
+        {
+          const std::array<EK::Point_3, 3>& t = new_triangles[ti].first;
+          visitor.new_subtriangle(offset+ti, tri_inter_ids_inverse[new_triangles[ti].second]);
+          triangle_buffer[ti] = CGAL::make_array(concurrent_get_point_id(t[0]), concurrent_get_point_id(t[1]), concurrent_get_point_id(t[2]));
+        }
+      }
+    );
+
+    // the map is now filled we can safely set the point ids
+    std::size_t pid_offset=soup_points.size();
+    soup_points.resize(pid_offset+iterators.size());
+#if ! defined(CGAL_NDEBUG) || defined(CGAL_DEBUG_PMP_AUTOREFINE)
+    exact_soup_points.resize(soup_points.size());
+#endif
+
+    tbb::parallel_for(tbb::blocked_range<size_t>(0, iterators.size()),
+      [&](const tbb::blocked_range<size_t>& r) {
+        for (size_t ti = r.begin(); ti != r.end(); ++ti)
+        {
+          soup_points[pid_offset+ti] = to_input(iterators[ti]->first);
+#if ! defined(CGAL_NDEBUG) || defined(CGAL_DEBUG_PMP_AUTOREFINE)
+          exact_soup_points[pid_offset+ti] = iterators[ti]->first;
+#endif
+          iterators[ti]->second=pid_offset+ti;
+        }
+      }
+    );
+
+    tbb::parallel_for(tbb::blocked_range<size_t>(0, new_triangles.size()),
+      [&](const tbb::blocked_range<size_t>& r) {
+          for (size_t ti = r.begin(); ti != r.end(); ++ti)
+          {
+            soup_triangles_out[offset + ti] =
+              { triangle_buffer[ti][0]->second,
+                triangle_buffer[ti][1]->second,
+                triangle_buffer[ti][2]->second };
+          }
+        }
+    );
+#endif
+  }
+  else
+#endif
+  {
+#ifdef CGAL_AUTOREF_USE_DEBUG_PARALLEL_TIMERS
+    mode = "sequential";
+#endif
+    soup_triangles_out.reserve(offset + new_triangles.size());
+    for (const std::pair<std::array<EK::Point_3,3>, std::size_t>& t_and_id : new_triangles)
+    {
+      visitor.new_subtriangle(soup_triangles_out.size(), tri_inter_ids_inverse[t_and_id.second]);
+      soup_triangles_out.push_back({ get_point_id(t_and_id.first[0]),
+                                     get_point_id(t_and_id.first[1]),
+                                     get_point_id(t_and_id.first[2]) });
+    }
+  }
+
+
+
+#ifdef CGAL_AUTOREF_USE_DEBUG_PARALLEL_TIMERS
+  t.stop();
+  std::cout << t.time() << " sec. for #4 (" << mode << ")" << std::endl;
+  t.reset();
+#endif
+
+#ifndef CGAL_NDEBUG
+  CGAL_PMP_AUTOREFINE_VERBOSE("check soup");
+  CGAL_assertion( !does_triangle_soup_self_intersect<Concurrency_tag>(exact_soup_points, soup_triangles_out) );
+#else
+#ifdef CGAL_DEBUG_PMP_AUTOREFINE
+  CGAL_PMP_AUTOREFINE_VERBOSE("check soup");
+  if (does_triangle_soup_self_intersect<Concurrency_tag>(exact_soup_points, soup_triangles_out))
+    throw std::runtime_error("ERROR: invalid output, there is most probably a bug");
+#endif
+#endif
+  using std::swap;
+  swap(soup_triangles, soup_triangles_out);
+
+  CGAL_PMP_AUTOREFINE_VERBOSE("done");
+}
+
+/**
+ * \ingroup PMP_corefinement_grp
+ * refines a triangle mesh so that no triangles intersects in their interior.
+ *
+ * Note that this function is only provided as a shortcut for calling `autorefine_triangle_soup()`
+ * with a mesh. For any advance usage the aforementioned function should be called directly.
+ *
+ * @tparam TriangleMesh a model of `HalfedgeListGraph`, `FaceListGraph`, and `MutableFaceGraph`
+ * @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
+ *
+ * @param tm input triangulated surface mesh
+ * @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+ *
+ * @warning `clear(tm)` will be called before filling `tm` with the refined mesh.
+ *
+ * \cgalNamedParamsBegin
+ *   \cgalParamNBegin{concurrency_tag}
+ *     \cgalParamDescription{a tag indicating if the task should be done using one or several threads.}
+ *     \cgalParamType{Either `CGAL::Sequential_tag`, or `CGAL::Parallel_tag`, or `CGAL::Parallel_if_available_tag`}
+ *     \cgalParamDefault{`CGAL::Sequential_tag`}
+ *   \cgalParamNEnd
+ *   \cgalParamNBegin{geom_traits}
+ *     \cgalParamDescription{an instance of a geometric traits class}
+ *     \cgalParamType{a class model of `PMPSelfIntersectionTraits`}
+ *     \cgalParamDefault{a \cgal Kernel deduced from the point type, using `CGAL::Kernel_traits`}
+ *     \cgalParamExtra{The geometric traits class must be compatible with the vertex point type.}
+ *   \cgalParamNEnd
+ *   \cgalParamNBegin{vertex_point_map}
+ *     \cgalParamDescription{a property map associating points to the vertices of `tm`}
+ *     \cgalParamType{a class model of `ReadablePropertyMap` with `boost::graph_traits<TriangleMesh>::%vertex_descriptor`
+ *                    as key type and `%Point_3` as value type}
+ *     \cgalParamDefault{`boost::get(CGAL::vertex_point, tm)`}
+ *     \cgalParamExtra{If this parameter is omitted, an internal property map for `CGAL::vertex_point_t`
+ *                     must be available in `TriangleMesh`.}
+ *   \cgalParamNEnd
+ * \cgalNamedParamsEnd
+ *
+ */
+template <class TriangleMesh,
+          class NamedParameters = parameters::Default_named_parameters>
+void
+autorefine(      TriangleMesh& tm,
+           const NamedParameters& np = parameters::default_values())
+{
+  using parameters::choose_parameter;
+  using parameters::get_parameter;
+
+  typedef typename GetGeomTraits<TriangleMesh, NamedParameters>::type GT;
+  // GT traits = choose_parameter<GT>(get_parameter(np, internal_np::geom_traits));
+
+  std::vector<typename GT::Point_3> soup_points;
+  std::vector<std::array<std::size_t, 3> > soup_triangles;
+
+  polygon_mesh_to_polygon_soup(tm, soup_points, soup_triangles, np);
+
+  autorefine_triangle_soup(soup_points, soup_triangles, np);
+
+  clear(tm);
+  repair_polygon_soup(soup_points, soup_triangles);
+
+  duplicate_non_manifold_edges_in_polygon_soup(soup_points, soup_triangles);
+  polygon_soup_to_polygon_mesh(soup_points, soup_triangles, tm);
+}
+
+
+} } // end of CGAL::Polygon_mesh_processing
+
+#ifdef CGAL_LINKED_WITH_TBB
+#ifdef CGAL_HAS_DEFINED_TBB_PREVIEW_CONCURRENT_ORDERED_CONTAINERS
+#undef TBB_PREVIEW_CONCURRENT_ORDERED_CONTAINERS
+#endif
+#endif
+
+#endif // CGAL_POLYGON_MESH_PROCESSING_AUTOREFINEMENT_H
diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt
index 222f04417c9..44e32386664 100644
--- a/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt
+++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt
@@ -106,6 +106,7 @@ if(TARGET CGAL::TBB_support)
   target_link_libraries(test_pmp_distance PUBLIC CGAL::TBB_support)
   target_link_libraries(orient_polygon_soup_test PUBLIC CGAL::TBB_support)
   target_link_libraries(self_intersection_surface_mesh_test PUBLIC CGAL::TBB_support)
+  target_link_libraries(test_autorefinement PUBLIC CGAL::TBB_support)
 else()
   message(STATUS "NOTICE: Intel TBB was not found. Tests will use sequential code.")
 endif()
diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_autorefinement.cmd b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_autorefinement.cmd
index a2ffe2c2539..8c848a2f106 100644
--- a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_autorefinement.cmd
+++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_autorefinement.cmd
@@ -1,28 +1,28 @@
-data-autoref/test_01.off 0 0 4 1 4 0
-data-autoref/test_02.off 1 5 13 1 5 0
-data-autoref/test_03.off 4 8 18 0 18 0
-data-autoref/test_04.off 1 5 17 0 17 0
-data-autoref/test_05.off 1 5 17 0 17 0
-data-autoref/test_06.off 3 55 141 1 76 0
-data-autoref/test_07.off 1 4 10 1 4 0
-data-autoref/test_08.off 1 31 87 1 52 0
-data-autoref/test_09.off 1 4 5 1 4 0
-data-autoref/test_10.off 1 3 13 0 13 0
-data-autoref/test_11.off 1 3 12 0 12 0
-data-autoref/test_12.off 2 2 11 1 11 0
-data-autoref/test_13.off 1 5 16 1 8 0
-data-autoref/test_14.off 1 5 16 1 12 0
-data-autoref/test_15.off 3 8 16 1 12 0
-data-autoref/test_16.off 1 2 6 1 4 0
-data-autoref/test_17.off 1 2 6 1 4 0
-data-autoref/triple_inter_exception/triple.off 0 0 0 0 0 1
-data-autoref/triple_inter_exception/cubes_cpln_1.off 0 0 0 0 0 1
-data-autoref/triple_inter_exception/cubes_cpln_2.off 0 0 0 0 0 1
-data-autoref/triple_inter_exception/cubes_cpln_3.off 0 0 0 0 0 1
-data-autoref/open_01.off 1 65 1377 1 1313 0
-data-autoref/open_02.off 1 33 595 1 562 0
-data-autoref/cpln_01.off 18 42 48 1 30 0
-data-autoref/cpln_02.off 28 56 40 1 24 0
-data-autoref/cpln_03.off 15 35 42 1 27 0
-data-autoref/four_cubes.off 12 94 184 1 78 0
-data-autoref/spiral.off 7 14 26 0 26 0
+data-autoref/test_01.off 0 0 4 1 4 0 4 2
+data-autoref/test_02.off 1 5 13 1 5 0 10 17
+data-autoref/test_03.off 4 8 18 0 18 0 14 20
+data-autoref/test_04.off 1 5 17 0 17 0 13 20
+data-autoref/test_05.off 1 5 17 0 17 0 13 20
+data-autoref/test_06.off 3 55 141 1 76 0 92 248
+data-autoref/test_07.off 1 4 10 1 4 0 8 12
+data-autoref/test_08.off 1 31 87 1 52 0 57 148
+data-autoref/test_09.off 1 4 5 1 4 0 4 4
+data-autoref/test_10.off 1 3 13 0 13 0 10 13
+data-autoref/test_11.off 1 3 12 0 12 0 9 12
+data-autoref/test_12.off 2 2 11 1 11 0 9 6
+data-autoref/test_13.off 1 5 16 1 8 0 16 22
+data-autoref/test_14.off 1 5 16 1 12 0 16 22
+data-autoref/test_15.off 3 8 16 1 12 0 16 24
+data-autoref/test_16.off 1 2 6 1 4 0 6 2
+data-autoref/test_17.off 1 2 6 1 4 0 6 2
+data-autoref/triple_inter_exception/triple.off 0 0 0 0 0 1 15 18
+data-autoref/triple_inter_exception/cubes_cpln_1.off 0 0 0 0 0 1 66 224
+data-autoref/triple_inter_exception/cubes_cpln_2.off 0 0 0 0 0 1 54 196
+data-autoref/triple_inter_exception/cubes_cpln_3.off 0 0 0 0 0 1 61 204
+data-autoref/open_01.off 1 65 1377 1 1313 0 1317 2622
+data-autoref/open_02.off 1 33 595 1 562 0 565 1056
+data-autoref/cpln_01.off 18 42 48 1 30 0 30 88
+data-autoref/cpln_02.off 28 56 40 1 24 0 24 72
+data-autoref/cpln_03.off 15 35 42 1 27 0 27 76
+data-autoref/four_cubes.off 12 94 184 1 78 0 106 352
+data-autoref/spiral.off 7 14 26 0 26 0 19 31
diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_autorefinement.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_autorefinement.cpp
index 9c908d7ec24..3af3402a5f2 100644
--- a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_autorefinement.cpp
+++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_autorefinement.cpp
@@ -3,6 +3,10 @@
 
 #include <CGAL/Polygon_mesh_processing/intersection.h>
 #include <CGAL/Polygon_mesh_processing/corefinement.h>
+#include <CGAL/Polygon_mesh_processing/autorefinement.h>
+#include <CGAL/Polygon_mesh_processing/self_intersections.h>
+
+#include <CGAL/IO/polygon_soup_io.h>
 
 #include <fstream>
 #include <sstream>
@@ -13,23 +17,64 @@ typedef CGAL::Surface_mesh<K::Point_3>             Mesh;
 namespace PMP = CGAL::Polygon_mesh_processing;
 
 template <class TriangleMesh>
-struct My_visitor :
+struct My_exp_visitor :
   public CGAL::Polygon_mesh_processing::Corefinement::Default_visitor<TriangleMesh>
 {
   void after_subface_creations(TriangleMesh&){++(*i);}
 
-  My_visitor()
+  My_exp_visitor()
     : i (new int(0) )
   {}
 
   std::shared_ptr<int> i;
 };
 
-void test(const char* fname, std::size_t nb_polylines, std::size_t total_nb_points,
-          std::size_t nb_vertices_after_autorefine, bool all_fixed, std::size_t nb_vertices_after_fix,
-          bool triple_intersection)
+struct My_visitor
 {
-  std::cout << "Running tests on " << fname << "\n";
+  My_visitor(std::size_t nb_input, std::size_t expected_nb_output)
+    : nb_input(nb_input)
+    , expected_nb_output(expected_nb_output)
+  {}
+
+  ~My_visitor()
+  {
+    for(std::size_t i=0; i<tgt_check.size(); ++i)
+    {
+      assert( tgt_check[i]==1 );
+    }
+  }
+  void number_of_output_triangles(std::size_t nbt)
+  {
+    tgt_check.assign(expected_nb_output, 0);
+    assert(nbt==expected_nb_output);
+  }
+
+  void verbatim_triangle_copy(std::size_t tgt_id, std::size_t src_id)
+  {
+    assert(src_id<nb_input);
+    assert(tgt_id<expected_nb_output);
+    assert(tgt_check.size()==expected_nb_output && tgt_check[tgt_id]==0);
+    tgt_check[tgt_id]=1;
+  }
+
+  void new_subtriangle(std::size_t tgt_id, std::size_t src_id)
+  {
+    assert(src_id<nb_input);
+    assert(tgt_id<expected_nb_output);
+    assert(tgt_check.size()==expected_nb_output && tgt_check[tgt_id]==0);
+    tgt_check[tgt_id]=1;
+  }
+
+  std::size_t nb_input;
+  std::size_t expected_nb_output;
+  std::vector<int> tgt_check;
+};
+
+void test_coref_based(const char* fname, std::size_t nb_polylines, std::size_t total_nb_points,
+                      std::size_t nb_vertices_after_autorefine, bool all_fixed, std::size_t nb_vertices_after_fix,
+                      bool triple_intersection)
+{
+  std::cout << "Running tests (coref based) on " << fname << "\n";
   std::ifstream input(fname);
 
   Mesh mesh;
@@ -41,7 +86,7 @@ void test(const char* fname, std::size_t nb_polylines, std::size_t total_nb_poin
   input.close();
   std::size_t nb_vertices_before_autorefine = num_vertices(mesh);
 
-// Testing surface_self_intersection()
+// Testing PMP::experimental::surface_self_intersection()
   try{
     std::vector< std::vector<K::Point_3> >polylines;
     PMP::experimental::surface_self_intersection(mesh, std::back_inserter(polylines));
@@ -57,9 +102,9 @@ void test(const char* fname, std::size_t nb_polylines, std::size_t total_nb_poin
     assert( triple_intersection );
   }
 
-// Testing autorefine()
+// Testing PMP::experimental::autorefine()
   try{
-    My_visitor<Mesh> visitor;
+    My_exp_visitor<Mesh> visitor;
     PMP::experimental::autorefine(mesh,
       CGAL::parameters::visitor(visitor));
     mesh.collect_garbage();
@@ -72,7 +117,7 @@ void test(const char* fname, std::size_t nb_polylines, std::size_t total_nb_poin
     assert( triple_intersection );
   }
 
-// Testing autorefine_and_remove_self_intersections()
+// Testing PMP::experimental::autorefine_and_remove_self_intersections()
   try{
     input.open(fname);
     mesh.clear();
@@ -89,10 +134,42 @@ void test(const char* fname, std::size_t nb_polylines, std::size_t total_nb_poin
   }
 }
 
+template <class Tag>
+void test(const char* fname, std::size_t nb_vertices_after_autorefine, std::size_t expected_nb_output, Tag tag)
+{
+  std::cout << "Running tests on " << fname;
+  if (std::is_same_v<Tag, CGAL::Sequential_tag>)
+    std::cout << " (Sequential)\n";
+  else
+    std::cout << " (Parallel)\n";
+
+  std::vector<K::Point_3> points;
+  std::vector< std::vector<std::size_t> > triangles;
+  if (!CGAL::IO::read_polygon_soup(fname, points, triangles))
+  {
+    std::cerr << "  Input mesh is not a valid file." << std::endl;
+    exit(EXIT_FAILURE);
+  }
+
+// Testing autorefine()
+  My_visitor visitor(triangles.size(), expected_nb_output);
+  PMP::autorefine_triangle_soup(points, triangles, CGAL::parameters::visitor(visitor).concurrency_tag(tag));
+  assert( nb_vertices_after_autorefine==points.size());
+  assert( expected_nb_output==triangles.size());
+  assert( !PMP::does_triangle_soup_self_intersect(points, triangles) );
+//  CGAL::IO::write_polygon_soup("/tmp/debug.off", points, triangles);
+}
+
 int main(int argc, const char** argv)
 {
   // file nb_polylines total_nb_points nb_vertices_after_autorefine all_fixed nb_vertices_after_fix triple_intersection
-  for (int i=0;i<(argc-1)/7; ++i)
-    test(argv[1+7*i], atoi(argv[1+7*i+1]), atoi(argv[1+7*i+2]),
-         atoi(argv[1+7*i+3]), atoi(argv[1+7*i+4])==0?false:true, atoi(argv[1+7*i+5]), atoi(argv[1+7*i+6])==0?false:true);
+  for (int i=0;i<(argc-1)/9; ++i)
+  {
+    test_coref_based(argv[1+9*i], atoi(argv[1+9*i+1]), atoi(argv[1+9*i+2]),
+                     atoi(argv[1+9*i+3]), atoi(argv[1+9*i+4])==0?false:true, atoi(argv[1+9*i+5]), atoi(argv[1+9*i+6])==0?false:true);
+    test(argv[1+9*i], atoi(argv[1+9*i+7]), atoi(argv[1+9*i+8]), CGAL::Sequential_tag());
+#ifdef CGAL_LINKED_WITH_TBB
+    test(argv[1+9*i], atoi(argv[1+9*i+7]), atoi(argv[1+9*i+8]), CGAL::Parallel_tag());
+#endif
+  }
 }
diff --git a/Polyhedron/demo/Polyhedron/Plugins/PMP/CMakeLists.txt b/Polyhedron/demo/Polyhedron/Plugins/PMP/CMakeLists.txt
index 2651b3601e7..232e6f364d6 100644
--- a/Polyhedron/demo/Polyhedron/Plugins/PMP/CMakeLists.txt
+++ b/Polyhedron/demo/Polyhedron/Plugins/PMP/CMakeLists.txt
@@ -133,7 +133,11 @@ target_link_libraries(
 
 qt6_wrap_ui( repairUI_FILES RemoveNeedlesDialog.ui SelfSnapDialog.ui)
 polyhedron_demo_plugin(repair_polyhedron_plugin Repair_polyhedron_plugin  ${repairUI_FILES} KEYWORDS PMP)
-target_link_libraries(repair_polyhedron_plugin PUBLIC scene_points_with_normal_item scene_surface_mesh_item)
+target_link_libraries(repair_polyhedron_plugin PUBLIC scene_points_with_normal_item scene_surface_mesh_item scene_polygon_soup_item)
+if(TARGET CGAL::TBB_support)
+  target_link_libraries(repair_polyhedron_plugin PUBLIC CGAL::TBB_support)
+endif()
+
 
 if(TARGET CGAL::Eigen3_support)
   qt6_wrap_ui(isotropicRemeshingUI_FILES Isotropic_remeshing_dialog.ui)
diff --git a/Polyhedron/demo/Polyhedron/Plugins/PMP/Repair_polyhedron_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/PMP/Repair_polyhedron_plugin.cpp
index c4c68388d8c..b0ac1f85996 100644
--- a/Polyhedron/demo/Polyhedron/Plugins/PMP/Repair_polyhedron_plugin.cpp
+++ b/Polyhedron/demo/Polyhedron/Plugins/PMP/Repair_polyhedron_plugin.cpp
@@ -1,6 +1,7 @@
 #include <QtCore/qglobal.h>
 
 #include "Scene_surface_mesh_item.h"
+#include "Scene_polygon_soup_item.h"
 #include "Scene_points_with_normal_item.h"
 #include <CGAL/Three/Scene_interface.h>
 #include <CGAL/Three/Polyhedron_demo_plugin_interface.h>
@@ -19,6 +20,8 @@
 #include <CGAL/Polygon_mesh_processing/repair_degeneracies.h>
 #include <CGAL/Polygon_mesh_processing/merge_border_vertices.h>
 #include <CGAL/Polygon_mesh_processing/internal/Snapping/snap.h>
+#include <CGAL/Polygon_mesh_processing/autorefinement.h>
+#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
 
 #include "ui_RemoveNeedlesDialog.h"
 #include "ui_SelfSnapDialog.h"
@@ -52,8 +55,9 @@ public:
     actionDuplicateNMVertices = new QAction(tr("Duplicate Non-Manifold Vertices"), mw);
     actionExtractNMVertices = new QAction(tr("Extract Non-Manifold Vertices"), mw);
     actionMergeDuplicatedVerticesOnBoundaryCycles = new QAction(tr("Merge Duplicated Vertices on Boundary Cycles"), mw);
-    actionAutorefine = new QAction(tr("Autorefine Mesh"), mw);
-    actionAutorefineAndRMSelfIntersections = new QAction(tr("Autorefine and Remove Self-Intersections"), mw);
+    actionAutorefine = new QAction(tr("Autorefine Mesh (Deprecated)"), mw);
+    actionNewAutorefine = new QAction(tr("Autorefine"), mw);
+    actionAutorefineAndRMSelfIntersections = new QAction(tr("Autorefine and Remove Self-Intersections (Deprecated)"), mw);
     actionRemoveNeedlesAndCaps = new QAction(tr("Remove Needles And Caps"));
     actionSnapBorders = new QAction(tr("Snap Boundaries"));
 
@@ -65,6 +69,7 @@ public:
     actionExtractNMVertices->setObjectName("actionExtractNMVertices");
     actionMergeDuplicatedVerticesOnBoundaryCycles->setObjectName("actionMergeDuplicatedVerticesOnBoundaryCycles");
     actionAutorefine->setObjectName("actionAutorefine");
+    actionNewAutorefine->setObjectName("actionNewAutorefine");
     actionAutorefineAndRMSelfIntersections->setObjectName("actionAutorefineAndRMSelfIntersections");
     actionRemoveNeedlesAndCaps->setObjectName("actionRemoveNeedlesAndCaps");
     actionSnapBorders->setObjectName("actionSnapBorders");
@@ -77,6 +82,7 @@ public:
     actionExtractNMVertices->setProperty("subMenuName", "Polygon Mesh Processing/Repair");
     actionMergeDuplicatedVerticesOnBoundaryCycles->setProperty("subMenuName", "Polygon Mesh Processing/Repair");
     actionAutorefine->setProperty("subMenuName", "Polygon Mesh Processing/Repair/Experimental");
+    actionNewAutorefine->setProperty("subMenuName", "Polygon Mesh Processing/Repair");
     actionAutorefineAndRMSelfIntersections->setProperty("subMenuName", "Polygon Mesh Processing/Repair/Experimental");
     actionSnapBorders->setProperty("subMenuName", "Polygon Mesh Processing/Repair/Experimental");
 
@@ -93,16 +99,29 @@ public:
                              << actionExtractNMVertices
                              << actionMergeDuplicatedVerticesOnBoundaryCycles
                              << actionAutorefine
+                             << actionNewAutorefine
                              << actionAutorefineAndRMSelfIntersections
                              << actionRemoveNeedlesAndCaps
                              << actionSnapBorders;
   }
 
-  bool applicable(QAction*) const
+  bool applicable(QAction* action) const
   {
-    int item_id = scene->mainSelectionIndex();
-    return qobject_cast<Scene_surface_mesh_item*>(scene->item(item_id));
+    if (action!=actionNewAutorefine)
+    {
+      int item_id = scene->mainSelectionIndex();
+      return qobject_cast<Scene_surface_mesh_item*>(scene->item(item_id));
+    }
+    for (Scene_interface::Item_id index : scene->selectionIndices())
+    {
+      if (qobject_cast<Scene_surface_mesh_item*>(scene->item(index)))
+        return true;
+      if (qobject_cast<Scene_polygon_soup_item*>(scene->item(index)))
+        return true;
+    }
+    return false;
   }
+
   template <typename Item>
   void on_actionRemoveIsolatedVertices_triggered(Scene_interface::Item_id index);
   template <typename Item>
@@ -120,6 +139,8 @@ public:
   template <typename Item>
   void on_actionAutorefine_triggered(Scene_interface::Item_id index);
   template <typename Item>
+  void on_actionNewAutorefine_triggered(const std::vector<Scene_interface::Item_id>& indices);
+  template <typename Item>
   void on_actionAutorefineAndRMSelfIntersections_triggered(Scene_interface::Item_id index);
 
 public Q_SLOTS:
@@ -131,6 +152,7 @@ public Q_SLOTS:
   void on_actionExtractNMVertices_triggered();
   void on_actionMergeDuplicatedVerticesOnBoundaryCycles_triggered();
   void on_actionAutorefine_triggered();
+  void on_actionNewAutorefine_triggered();
   void on_actionAutorefineAndRMSelfIntersections_triggered();
   void on_actionRemoveNeedlesAndCaps_triggered();
   void on_actionSnapBorders_triggered();
@@ -144,6 +166,7 @@ private:
   QAction* actionExtractNMVertices;
   QAction* actionMergeDuplicatedVerticesOnBoundaryCycles;
   QAction* actionAutorefine;
+  QAction* actionNewAutorefine;
   QAction* actionAutorefineAndRMSelfIntersections;
   QAction* actionRemoveNeedlesAndCaps;
   QAction* actionSnapBorders;
@@ -421,6 +444,78 @@ void Polyhedron_demo_repair_polyhedron_plugin::on_actionAutorefine_triggered()
   QApplication::restoreOverrideCursor();
 }
 
+template <typename Item>
+void Polyhedron_demo_repair_polyhedron_plugin::on_actionNewAutorefine_triggered(const std::vector<Scene_interface::Item_id>& indices)
+{
+  namespace PMP = CGAL::Polygon_mesh_processing;
+  Polygon_soup::Points points;
+  Polygon_soup::Polygons polygons;
+
+  if (indices.size()==1)
+  {
+    if (Scene_surface_mesh_item* smi_ptr = qobject_cast<Scene_surface_mesh_item*>(scene->item(indices[0])))
+      PMP::polygon_mesh_to_polygon_soup(*smi_ptr->polyhedron(), points, polygons);
+    else if (Scene_polygon_soup_item* spi_ptr = qobject_cast<Scene_polygon_soup_item*>(scene->item(indices[0])))
+    {
+      points = spi_ptr->points();
+      polygons = spi_ptr->polygons();
+    }
+  }
+  else
+  {
+    for (Scene_interface::Item_id id : indices)
+    {
+      Polygon_soup::Points l_points;
+      Polygon_soup::Polygons l_polygons;
+
+      if (Scene_surface_mesh_item* smi_ptr = qobject_cast<Scene_surface_mesh_item*>(scene->item(id)))
+        PMP::polygon_mesh_to_polygon_soup(*smi_ptr->polyhedron(), l_points, l_polygons);
+      else if (Scene_polygon_soup_item* spi_ptr = qobject_cast<Scene_polygon_soup_item*>(scene->item(id)))
+      {
+        l_points = spi_ptr->points();
+        l_polygons = spi_ptr->polygons();
+      }
+      std::size_t offset=points.size();
+      points.insert(points.end(), l_points.begin(), l_points.end());
+      std::size_t psize=polygons.size();
+      polygons.insert(polygons.end(), l_polygons.begin(), l_polygons.end());
+      for (std::size_t i=psize; i<polygons.size(); ++i)
+        for(std::size_t& id : polygons[i])
+          id+=offset;
+    }
+  }
+
+  PMP::triangulate_polygons(points, polygons);
+  PMP::autorefine_triangle_soup(points, polygons,
+                                CGAL::parameters::concurrency_tag(CGAL::Parallel_if_available_tag()));
+
+  Scene_polygon_soup_item* new_item = new Scene_polygon_soup_item();
+  new_item->load(points, polygons);
+  QString name = scene->item(indices[0])->name();
+  for (std::size_t k=1; k<indices.size(); ++k)
+    name += " + " + scene->item(indices[k])->name();
+  new_item->setName(name+" autorefined");
+
+  scene->addItem(new_item);
+  new_item->invalidateOpenGLBuffers();
+  Q_EMIT new_item->itemChanged();
+}
+
+void Polyhedron_demo_repair_polyhedron_plugin::on_actionNewAutorefine_triggered()
+{
+  std::vector<Scene_interface::Item_id> indices;
+  for (Scene_interface::Item_id index : scene->selectionIndices())
+  {
+    if (qobject_cast<Scene_surface_mesh_item*>(scene->item(index)))
+      indices.push_back(index);
+    else if (qobject_cast<Scene_polygon_soup_item*>(scene->item(index)))
+      indices.push_back(index);
+  }
+  QApplication::setOverrideCursor(Qt::WaitCursor);
+  on_actionNewAutorefine_triggered<Scene_surface_mesh_item>(indices);
+  QApplication::restoreOverrideCursor();
+}
+
 template <typename Item>
 void Polyhedron_demo_repair_polyhedron_plugin::on_actionAutorefineAndRMSelfIntersections_triggered(Scene_interface::Item_id index)
 {
diff --git a/Property_map/include/CGAL/property_map.h b/Property_map/include/CGAL/property_map.h
index 735c1f5da49..1154ab87b42 100644
--- a/Property_map/include/CGAL/property_map.h
+++ b/Property_map/include/CGAL/property_map.h
@@ -661,7 +661,7 @@ struct Boolean_property_map
     return pm.set_ptr->count(k) != 0;
   }
 
-  friend void put(Boolean_property_map<Set>& pm, const key_type& k, bool v)
+  friend void put(Boolean_property_map<Set> pm, const key_type& k, bool v)
   {
     CGAL_assertion(pm.set_ptr!=nullptr);
     if (v)
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 51e1856e754..9b172f5b4c1 100644
--- a/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h
+++ b/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h
@@ -47,6 +47,7 @@ CGAL_add_named_parameter(implementation_tag_t, implementation_tag, implementatio
 CGAL_add_named_parameter(prevent_unselection_t, prevent_unselection, prevent_unselection)
 
 CGAL_add_named_parameter(verbose_t, verbose, verbose)
+CGAL_add_named_parameter(concurrency_tag_t, concurrency_tag, concurrency_tag)
 
 // List of named parameters used for IO
 CGAL_add_named_parameter(vertex_normal_output_iterator_t, vertex_normal_output_iterator, vertex_normal_output_iterator)
diff --git a/Triangulation_2/include/CGAL/Constrained_triangulation_2.h b/Triangulation_2/include/CGAL/Constrained_triangulation_2.h
index 048b373187e..3b00bd22e5f 100644
--- a/Triangulation_2/include/CGAL/Constrained_triangulation_2.h
+++ b/Triangulation_2/include/CGAL/Constrained_triangulation_2.h
@@ -330,11 +330,11 @@ public:
 #if 1
   template <class Segment_2>
   static decltype(auto) get_source(const Segment_2& segment){
-    return segment.source();
+    return segment[0];
   }
   template <class Segment_2>
   static decltype(auto) get_target(const Segment_2& segment){
-    return segment.target();
+    return segment[1];
   }
 
   static const Point& get_source(const Constraint& cst){