diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/orient_polygon_soup.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/orient_polygon_soup.h
index 6a17c34bce7..f815cf2d8bb 100644
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/orient_polygon_soup.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/orient_polygon_soup.h
@@ -21,8 +21,17 @@
 #include <CGAL/tuple.h>
 #include <CGAL/array.h>
 #include <CGAL/assertions.h>
+#include <CGAL/AABB_tree.h>
+#include <CGAL/AABB_face_graph_triangle_primitive.h>
+#include <CGAL/AABB_traits.h>
+
+#include <CGAL/Polygon_mesh_processing/shape_predicates.h>
+#include <CGAL/Polygon_mesh_processing/compute_normal.h>
+
 #include <boost/container/flat_set.hpp>
 #include <boost/container/flat_map.hpp>
+#include <boost/iterator/counting_iterator.hpp>
+#include <boost/iterator/filter_iterator.hpp>
 
 #include <set>
 #include <map>
@@ -482,6 +491,122 @@ bool orient_polygon_soup(PointRange& points,
   return inital_nb_pts==points.size();
 }
 
+/*!
+ * Duplicate each point \a p at which the intersection
+ * of an infinitesimally small ball centered at \a p
+ * with the polygons incident to it is not a topological disk.
+ *
+ * @tparam PointRange a model of the concepts `RandomAccessContainer`
+ * and `BackInsertionSequence` whose value type is the point type.
+ * @tparam PolygonRange a model of the concept `RandomAccessContainer`
+ * whose value_type is a model of the concept `RandomAccessContainer`
+ * whose value_type is `std::size_t`.
+ *
+ * @param points points of the soup of polygons. Some additional points might be pushed back to resolve
+ *               non-manifoldness or non-orientability issues.
+ * @param polygons each element in the vector describes a polygon using the index of the points in `points`.
+ *                 If needed the order of the indices of a polygon might be reversed.
+ * @return `true`  if the orientation operation succeded.
+ * @return `false` if some points were duplicated, thus producing a self-intersecting polyhedron.
+ * @sa `orient_polygon_soup()`
+ */
+template <class PointRange, class PolygonRange>
+bool
+duplicate_incompatible_edges_in_polygon_soup(PointRange& points,
+                                            PolygonRange& polygons)
+{
+  std::size_t inital_nb_pts = points.size();
+  typedef CGAL::Polygon_mesh_processing::internal::
+    Polygon_soup_orienter<PointRange, PolygonRange> Orienter;
+
+  Orienter orienter(points, polygons);
+  orienter.fill_edge_map();
+  // make edges to duplicate
+  for(std::size_t i1=0;i1<points.size();++i1)
+    for(const typename Orienter::Internal_map_type::value_type& i2_and_pids : orienter.edges[i1])
+      if (i2_and_pids.second.size() > 1)
+        orienter.set_edge_marked(i1,i2_and_pids.first,orienter.marked_edges);
+  orienter.duplicate_singular_vertices();
+
+  return inital_nb_pts==points.size();
+}
+
+/*!
+ * Orient each triangle of a triangle soup using the orientation of its
+ * closest non degenerate triangle in `tm_ref`.
+ * \tparam Concurrency_tag enables sequential versus parallel orientation.
+                        Possible values are `Sequential_tag` (the default) and
+                        `Parallel_tag`.
+ * \tparam Point_3 the point type of the soup.
+ * @tparam TriangleRange a model of the concept `RandomAccessContainer`
+ * whose value_type is a model of the concept `RandomAccessContainer`
+ * whose value_type is `std::size_t`and of size 3.
+ * @tparam TriangleMesh a model of `FaceListGraph` and `MutableFaceGraph` .
+ *
+ * \param tm_ref the reference TriangleMesh.
+ * \param points the points of the soup.
+ * \param triangles the triangles of the soup.
+ */
+template <class Concurrency_tag, class Point_3, class TriangleRange, class TriangleMesh>
+void
+orient_triangle_soup_with_reference_triangle_mesh(
+  const TriangleMesh& tm_ref,
+  std::vector<Point_3>& points,
+  TriangleRange& triangles)
+{
+  namespace PMP = CGAL::Polygon_mesh_processing;
+
+  typedef boost::graph_traits<TriangleMesh> GrT;
+  typedef typename GrT::face_descriptor face_descriptor;
+  typedef typename Kernel_traits<Point_3>::Kernel K;
+  typedef std::function<bool(face_descriptor)> Face_predicate;
+  Face_predicate is_not_deg =
+    [&tm_ref](face_descriptor f)
+    {
+      return !PMP::is_degenerate_triangle_face(f, tm_ref);
+    };
+
+  // build a tree filtering degenerate faces
+  typedef CGAL::AABB_face_graph_triangle_primitive<TriangleMesh> Primitive;
+  typedef CGAL::AABB_traits<K, Primitive> Tree_traits;
+
+  boost::filter_iterator<Face_predicate, typename GrT::face_iterator>
+    begin(is_not_deg, faces(tm_ref).begin(), faces(tm_ref).end()),
+    end(is_not_deg, faces(tm_ref).end(), faces(tm_ref).end());
+
+  CGAL::AABB_tree<Tree_traits> tree(begin, end, tm_ref);
+
+  // now orient the faces
+  tree.build();
+  tree.accelerate_distance_queries();
+  auto process_facet =
+    [&points, &tree, &tm_ref, &triangles](std::size_t fid) {
+      const auto& p0 = points[triangles[fid][0]];
+      const auto& p1 = points[triangles[fid][1]];
+      const auto& p2 = points[triangles[fid][2]];
+      const typename K::Point_3 mid = CGAL::centroid(p0, p1, p2);
+      std::pair<typename K::Point_3, face_descriptor> pt_and_f =
+        tree.closest_point_and_primitive(mid);
+      auto face_ref_normal = PMP::compute_face_normal(pt_and_f.second, tm_ref);
+      if(face_ref_normal * cross_product(p1-p0, p2-p0) < 0) {
+        std::swap(triangles[fid][1], triangles[fid][2]);
+      }
+    };
+
+#if !defined(CGAL_LINKED_WITH_TBB)
+  CGAL_static_assertion_msg (!(boost::is_convertible<Concurrency_tag, CGAL::Parallel_tag>::value),
+                             "Parallel_tag is enabled but TBB is unavailable.");
+#else
+  if (boost::is_convertible<Concurrency_tag,CGAL::Parallel_tag>::value)
+    tbb::parallel_for(std::size_t(0), triangles.size(), std::size_t(1), process_facet);
+  else
+#endif
+    std::for_each(
+      boost::counting_iterator<std::size_t> (0),
+      boost::counting_iterator<std::size_t> (triangles.size()),
+      process_facet);
+}
+
 } }//end namespace CGAL::Polygon_mesh_processing
 
 #include <CGAL/enable_warnings.h>
diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/orientation.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/orientation.h
index 182a9a866f7..e4d19bb268c 100644
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/orientation.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/orientation.h
@@ -19,6 +19,7 @@
 
 #include <algorithm>
 #include <CGAL/Polygon_mesh_processing/connected_components.h>
+#include <CGAL/Polygon_mesh_processing/stitch_borders.h>
 #include <CGAL/Polygon_mesh_processing/compute_normal.h>
 #include <CGAL/Polygon_mesh_processing/internal/named_function_params.h>
 #include <CGAL/Polygon_mesh_processing/internal/named_params_helper.h>
@@ -646,6 +647,119 @@ void orient_to_bound_a_volume(TriangleMesh& tm)
 {
   orient_to_bound_a_volume(tm, parameters::all_default());
 }
+
+/*!
+ * \param tm
+ * \todo add namedparameters and fidmap management
+ */
+template <class TriangleMesh>
+void merge_reversible_connected_components(TriangleMesh& tm)
+{
+
+  namespace PMP = CGAL::Polygon_mesh_processing;
+
+  typedef boost::graph_traits<TriangleMesh> GrT;
+  typedef typename GrT::face_descriptor face_descriptor;
+  typedef typename GrT::halfedge_descriptor halfedge_descriptor;
+
+  typedef typename boost::property_map<TriangleMesh, CGAL::vertex_point_t >::type Vpm;
+  typedef typename boost::property_traits<Vpm>::value_type Point_3;
+  Vpm vpm = get(CGAL::vertex_point, tm);
+
+  typedef typename boost::property_map<TriangleMesh, CGAL::dynamic_face_property_t<std::size_t> >::type Fidmap;
+  Fidmap fim = get(CGAL::dynamic_face_property_t<std::size_t>(), tm);
+  std::size_t i=0;
+  for (face_descriptor f : faces(tm))
+    put(fim, f, i++);
+
+  Fidmap f_cc_ids = get(CGAL::dynamic_face_property_t<std::size_t>(), tm);
+  std::size_t nb_cc = PMP::connected_components(tm, f_cc_ids, CGAL::parameters::face_index_map(fim));
+
+  std::vector<std::size_t> nb_faces_per_cc(nb_cc, 0);
+  for (face_descriptor f : faces(tm))
+    nb_faces_per_cc[ get(f_cc_ids, f) ]+=1;
+
+  std::map< std::pair<Point_3, Point_3>, std::vector<halfedge_descriptor> > border_hedges_map;
+  std::vector<halfedge_descriptor> border_hedges;
+  typedef typename boost::property_map<TriangleMesh, CGAL::dynamic_halfedge_property_t<std::size_t> >::type Hidmap;
+  Hidmap him = get(CGAL::dynamic_halfedge_property_t<std::size_t>(), tm);
+  const std::size_t DV(-1);
+
+  // fill endpoints -> hedges
+  for (halfedge_descriptor h : halfedges(tm))
+  {
+    if ( CGAL::is_border(h, tm) )
+    {
+      put(him, h, DV);
+      border_hedges_map[std::make_pair(get(vpm, source(h, tm)), get(vpm, target(h, tm)))].push_back(h);
+      border_hedges.push_back(h);
+    }
+  }
+
+  // set the id of boundary cc
+  std::size_t id=0;
+  for(halfedge_descriptor h : border_hedges)
+  {
+    if (get(him,h) == DV)
+    {
+      for (halfedge_descriptor hh : CGAL::halfedges_around_face(h, tm))
+      {
+        put(him, hh, id);
+      }
+      ++id;
+    }
+  }
+
+  // check boundary fully matching
+  std::vector<bool> border_cycle_handled(id, false);
+  std::set<std::size_t> ccs_to_reverse;
+  for ( const auto& p : border_hedges_map )
+  {
+    const std::vector<halfedge_descriptor>& hedges = p.second;
+    if ( hedges.size()==2 && !border_cycle_handled[ get(him, hedges[0]) ]
+                          && !border_cycle_handled[ get(him, hedges[1]) ] )
+    {
+      border_cycle_handled[ get(him, hedges[0]) ] = true;
+      border_cycle_handled[ get(him, hedges[1]) ] = true;
+
+      halfedge_descriptor h2=hedges[1];
+      bool did_break=false;
+      for(halfedge_descriptor h1 : CGAL::halfedges_around_face(hedges[0], tm))
+      {
+        if (get(vpm, target(h1,tm)) != get(vpm, target(h2,tm)))
+        {
+          did_break=true;
+          break;
+        }
+        h2=next(h2,tm);
+      }
+      if (!did_break && h2==hedges[1])
+      {
+        std::size_t cc_id1 = get(f_cc_ids, face(opposite(hedges[0], tm), tm)),
+                    cc_id2 = get(f_cc_ids, face(opposite(hedges[1], tm), tm));
+
+        if (cc_id1!=cc_id2)
+        {
+          if (ccs_to_reverse.count(cc_id1)==0 && ccs_to_reverse.count(cc_id2)==0)
+            ccs_to_reverse.insert( nb_faces_per_cc[cc_id1] < nb_faces_per_cc[cc_id2] ? cc_id1 : cc_id2 );
+        }
+      }
+    }
+  }
+
+  // reverse ccs and stitches boundaries
+  std::vector<face_descriptor> faces_to_reverse;
+  for (face_descriptor f : faces(tm))
+    if ( ccs_to_reverse.count( get(f_cc_ids, f) ) != 0 )
+      faces_to_reverse.push_back(f);
+
+  if ( !faces_to_reverse.empty() )
+  {
+    PMP::reverse_face_orientations(faces_to_reverse, tm);
+    PMP::stitch_borders(tm);
+  }
+}
+
 } // namespace Polygon_mesh_processing
 } // namespace CGAL
 #endif // CGAL_ORIENT_POLYGON_MESH_H