diff --git a/Installation/CHANGES.md b/Installation/CHANGES.md
index cfc731b036f..3badada23f3 100644
--- a/Installation/CHANGES.md
+++ b/Installation/CHANGES.md
@@ -51,9 +51,8 @@ Release date: June 2022
 
 ### [Polygon Mesh Processing](https://doc.cgal.org/5.5/Manual/packages.html#PkgPolygonMeshProcessing)
 -   Added the function `CGAL::Polygon_mesh_processing::orient_triangle_soup_with_reference_triangle_soup()`, which enables re-orienting the faces of a triangle soup based on the orientation of the nearest face in a reference triangle soup.
--   Added the function `CGAL::Polygon_mesh_processing::tangential_relaxation()`, which
-applies an area-based tangential mesh smoothing to the vertices of a surface triangle mesh.
-
+-   Added the function `CGAL::Polygon_mesh_processing::compatible_orientations()`, which enables to retrieve the (in)compatibility of orientations of faces from different connected components.
+-   Added the function `CGAL::Polygon_mesh_processing::tangential_relaxation()`, which applies an area-based tangential mesh smoothing to the vertices of a surface triangle mesh.
 
 ### [2D Polygons](https://doc.cgal.org/5.5/Manual/packages.html#PkgPolygon2)
 
diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt
index 014ccb58919..6281ae0cdba 100644
--- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt
+++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt
@@ -132,6 +132,7 @@ The page \ref bgl_namedparameters "Named Parameters" describes their usage.
 - `CGAL::Polygon_mesh_processing::orient_triangle_soup_with_reference_triangle_mesh()`
 - `CGAL::Polygon_mesh_processing::orient_triangle_soup_with_reference_triangle_soup()`
 - `CGAL::Polygon_mesh_processing::merge_reversible_connected_components()`
+- `CGAL::Polygon_mesh_processing::compatible_orientations()`
 
 \cgalCRPSection{Hole Filling Functions}
 - `CGAL::Polygon_mesh_processing::triangulate_hole()`
diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt
index 1e172cdaa69..59cb8a647e0 100644
--- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt
+++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt
@@ -37,4 +37,5 @@
 \example Polygon_mesh_processing/orientation_pipeline_example.cpp
 \example Polygon_mesh_processing/triangulate_faces_split_visitor_example.cpp
 \example Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp
+\example Polygon_mesh_processing/cc_compatible_orientations.cpp
 */
diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt
index c83b338a588..04b0846cf9f 100644
--- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt
+++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt
@@ -94,6 +94,7 @@ create_single_source_cgal_program("orientation_pipeline_example.cpp")
 #create_single_source_cgal_program( "self_snapping_example.cpp")
 #create_single_source_cgal_program( "snapping_example.cpp")
 create_single_source_cgal_program("match_faces.cpp")
+create_single_source_cgal_program("cc_compatible_orientations.cpp")
 
 if(OpenMesh_FOUND)
 
diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/cc_compatible_orientations.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/cc_compatible_orientations.cpp
new file mode 100644
index 00000000000..8a47239c5aa
--- /dev/null
+++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/cc_compatible_orientations.cpp
@@ -0,0 +1,75 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Surface_mesh.h>
+
+#include <CGAL/Polygon_mesh_processing/orientation.h>
+#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
+#include <CGAL/Polygon_mesh_processing/stitch_borders.h>
+#include <CGAL/boost/graph/IO/polygon_mesh_io.h>
+
+#include <iostream>
+#include <string>
+#include <cstdlib>
+
+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;
+
+
+// create a mesh with many connected connected components that should
+// be reoriented to define a valid closed mesh
+void create_mesh_with_cc_to_orient(Mesh& mesh)
+{
+  const std::string filename = CGAL::data_file_path("meshes/elephant.off");
+  CGAL::IO::read_polygon_mesh(filename, mesh);
+
+  // turn the mesh into a triangle soup, duplicating all the vertices and shuffling orientations
+  std::vector<Point> points;
+  std::vector< std::array<std::size_t, 3> > triangles;
+  triangles.reserve(faces(mesh).size());
+  points.reserve(3*triangles.size());
+  for (Mesh::Face_index f : mesh.faces())
+  {
+    Mesh::Halfedge_index h = mesh.halfedge(f);
+    std::size_t s = points.size();
+    points.push_back(mesh.point(source(h,mesh)));
+    points.push_back(mesh.point(target(h,mesh)));
+    points.push_back(mesh.point(target(mesh.next(h),mesh)));
+    triangles.push_back( {s, s+1, s+2} );
+    if (std::rand() % 2 == 0)
+      std::swap(triangles.back()[0], triangles.back()[1]);
+  }
+
+  // load the soup into the mesh;
+  mesh.clear();
+  PMP::polygon_soup_to_polygon_mesh(points, triangles, mesh);
+}
+
+int main()
+{
+  Mesh mesh;
+  create_mesh_with_cc_to_orient(mesh);
+  CGAL::IO::write_polygon_mesh("to_orient.off", mesh, CGAL::parameters::stream_precision(17));
+
+  // determine face orientations to be reversed to create compatibility
+  auto fbm = mesh.add_property_map<Mesh::Face_index, bool>("fbm", false).first;
+  bool is_orientable = PMP::compatible_orientations(mesh, fbm);
+  assert(is_orientable);
+
+  // reverse orientation of faces with bit 1
+  std::vector<Mesh::Face_index> faces_to_reverse;
+  for (Mesh::Face_index f : mesh.faces())
+    if (get(fbm, f))
+      faces_to_reverse.push_back(f);
+  PMP::reverse_face_orientations(faces_to_reverse, mesh);
+
+  // there are still borders between previously incompatible faces: stitch to close the mesh
+  PMP::stitch_borders(mesh);
+
+  assert(CGAL::is_closed(mesh));
+  CGAL::IO::write_polygon_mesh("oriented_and_stitched.off", mesh, CGAL::parameters::stream_precision(17));
+
+  return 0;
+}
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 290c65ccebe..135d80edbcc 100644
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/orientation.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/orientation.h
@@ -1604,6 +1604,239 @@ void merge_reversible_connected_components(PolygonMesh& pm,
   }
 }
 
+
+/*!
+ * \ingroup PMP_orientation_grp
+ *
+ * identifies faces whose orientation must be reversed in order to enable stitching of connected components.
+ * Each face is assigned a bit (`false` or `true`)
+ * such that two faces have compatible orientations iff they are assigned the same bits.
+ *
+ * @tparam PolygonMesh a model of `HalfedgeListGraph`, `FaceGraph`.
+ * @tparam FaceBitMap a model of `WritablePropertyMap` with `face_descriptor` as key and `bool` as value_type
+ * @tparam NamedParameters a sequence of \ref bgl_namedparameters
+ *
+ * @param pm a surface mesh
+ * @param fbm face bit map indicating if a face orientation should be reversed to be stitchable
+ *        (see `CGAL::Polygon_mesh_processing::stitch_borders()`) with another face. If `false` is
+ *        returned, the map will not be filled.
+ * @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+ *
+ * @return `true` if `pm` can be reoriented and `false` otherwise.
+ *
+ * \cgalNamedParamsBegin
+ *   \cgalParamNBegin{vertex_point_map}
+ *     \cgalParamDescription{a property map associating points to the vertices of `pm`}
+ *     \cgalParamType{a class model of `ReadablePropertyMap` with `boost::graph_traits<PolygonMesh>::%vertex_descriptor`
+ *                    as key type and `%Point_3` as value type}
+ *     \cgalParamDefault{`boost::get(CGAL::vertex_point, pm)`}
+ *     \cgalParamExtra{If this parameter is omitted, an internal property map for `CGAL::vertex_point_t`
+ *                     should be available for the vertices of `pm`.}
+ *   \cgalParamNEnd
+ * \cgalNamedParamsEnd
+ *
+ * \sa reverse_face_orientations()
+ * \sa stitch_borders()
+ *
+ */
+template <class PolygonMesh, class FaceBitMap, class NamedParameters = parameters::Default_named_parameters>
+bool compatible_orientations(const PolygonMesh& pm,
+                             FaceBitMap fbm,
+                             const NamedParameters& np = parameters::default_values())
+{
+  typedef boost::graph_traits<PolygonMesh> GrT;
+  typedef typename GrT::face_descriptor face_descriptor;
+  typedef typename GrT::halfedge_descriptor halfedge_descriptor;
+
+  typedef typename GetVertexPointMap<PolygonMesh, NamedParameters>::const_type Vpm;
+
+  typedef typename boost::property_traits<Vpm>::value_type Point_3;
+  Vpm vpm = parameters::choose_parameter(parameters::get_parameter(np, internal_np::vertex_point),
+                                         get_const_property_map(vertex_point, pm));
+
+  typedef std::size_t F_cc_id; // Face cc-id
+  typedef std::size_t E_id; // Edge id
+
+  typedef dynamic_face_property_t<F_cc_id>                   Face_property_tag;
+  typedef typename boost::property_map<PolygonMesh, Face_property_tag>::const_type   Face_cc_map;
+  Face_cc_map f_cc_ids  = get(Face_property_tag(), pm);
+  F_cc_id nb_cc = connected_components(pm, f_cc_ids);
+
+  std::vector<std::size_t> nb_faces_per_cc(nb_cc, 0);
+  for (face_descriptor f : faces(pm))
+    nb_faces_per_cc[ get(f_cc_ids, f) ]+=1;
+
+  // collect border halfedges
+  std::vector<halfedge_descriptor> border_hedges;
+  for (halfedge_descriptor h : halfedges(pm))
+    if ( is_border(h, pm) )
+      border_hedges.push_back(h);
+  std::size_t nb_bh=border_hedges.size();
+
+  // compute the edge id of all border halfedges
+  typedef std::map< std::pair<Point_3, Point_3>, E_id> E_id_map;
+  E_id_map e_id_map;
+  E_id e_id = 0;
+
+  std::vector<E_id> eids;
+  eids.reserve(nb_bh);
+  for (halfedge_descriptor h : border_hedges)
+  {
+    std::pair< typename E_id_map::iterator, bool > insert_res =
+        e_id_map.insert(
+          std::make_pair(
+            make_sorted_pair(get(vpm, source(h, pm)),
+                             get(vpm, target(h,pm))), e_id) );
+    if (insert_res.second)
+      ++e_id;
+    eids.push_back(insert_res.first->second);
+  }
+
+  // fill incidence per edge
+  std::vector< std::vector<halfedge_descriptor> > incident_ccs_per_edge(e_id);
+  for (std::size_t i=0; i<nb_bh; ++i)
+    incident_ccs_per_edge[ eids[i] ].push_back(border_hedges[i]);
+
+  std::vector< std::vector<F_cc_id> > compatible_patches(nb_cc);
+  std::vector< std::vector<F_cc_id> > incompatible_patches(nb_cc);
+
+  for (std::vector<halfedge_descriptor>& v : incident_ccs_per_edge)
+  {
+    // ignore non-manifold edges
+    if (v.size()!=2) continue;
+    F_cc_id front_id=get(f_cc_ids, face(opposite(v.front(), pm), pm));
+    F_cc_id back_id=get(f_cc_ids, face(opposite(v.back(), pm), pm));
+
+    if (front_id==back_id) continue;
+
+    if (get(vpm, source(v.front(), pm))==get(vpm, target(v.back(), pm)))
+    {
+      compatible_patches[front_id].push_back(back_id);
+      compatible_patches[back_id].push_back(front_id);
+    }
+    else
+    {
+      incompatible_patches[front_id].push_back(back_id);
+      incompatible_patches[back_id].push_back(front_id);
+    }
+  }
+
+  for(F_cc_id cc_id=0; cc_id<nb_cc; ++cc_id)
+  {
+    std::sort(compatible_patches[cc_id].begin(), compatible_patches[cc_id].end());
+    std::sort(incompatible_patches[cc_id].begin(), incompatible_patches[cc_id].end());
+  }
+
+  // sort the connected components with potential matches using their number
+  // of faces (sorted by decreasing number of faces)
+  std::vector<bool> cc_bits(nb_cc, false);
+  std::vector<bool> cc_handled(nb_cc, false);
+
+  std::set< F_cc_id, std::function<bool(F_cc_id,F_cc_id)> > sorted_ids(
+    [&nb_faces_per_cc](F_cc_id i, F_cc_id j)
+    {return nb_faces_per_cc[i]==nb_faces_per_cc[j] ? i<j : nb_faces_per_cc[i]>nb_faces_per_cc[j];}
+  );
+  for(F_cc_id cc_id=0; cc_id<nb_cc; ++cc_id)
+    sorted_ids.insert(cc_id);
+
+  // consider largest CC first, default and set its bit to 0
+  for(F_cc_id cc_id : sorted_ids)
+  {
+    if (cc_handled[cc_id]) continue;
+
+    // extract compatible components
+    std::set<F_cc_id> bit_0_cc_set;
+    std::set<F_cc_id> bit_1_cc_set;
+    bit_0_cc_set.insert(cc_id);
+    std::vector<F_cc_id> stack_0=compatible_patches[cc_id];
+    std::vector<F_cc_id> stack_1=incompatible_patches[cc_id];
+
+    while( !stack_0.empty() || !stack_1.empty())
+    {
+      // increase the set of patches for bit 0 using compatible_patches
+      while( !stack_0.empty() )
+      {
+        F_cc_id back=stack_0.back();
+        stack_0.pop_back();
+        if (!bit_0_cc_set.insert(back).second) continue;
+        stack_0.insert(stack_0.end(), compatible_patches[back].begin(), compatible_patches[back].end());
+      }
+
+      // extract incompatible components
+      for (F_cc_id cid : bit_0_cc_set)
+        stack_1.insert(stack_1.end(), incompatible_patches[cid].begin(), incompatible_patches[cid].end());
+      // increase the set of patches for bit 1 using compatible_patches
+      while( !stack_1.empty() )
+      {
+        F_cc_id back=stack_1.back();
+        stack_1.pop_back();
+        if (!bit_1_cc_set.insert(back).second) continue;
+        stack_1.insert(stack_1.end(), compatible_patches[back].begin(), compatible_patches[back].end());
+      }
+      for (F_cc_id cid1 : bit_1_cc_set)
+        for (F_cc_id cid0 : incompatible_patches[cid1])
+          if( bit_0_cc_set.count(cid0)==0 )
+            stack_0.push_back(cid0);
+    }
+
+    // set intersection should be empty
+    std::vector<F_cc_id> inter;
+    std::set_intersection( bit_0_cc_set.begin(), bit_0_cc_set.end(),
+                           bit_1_cc_set.begin(), bit_1_cc_set.end(),
+                           std::back_inserter(inter));
+    if (!inter.empty())
+    {
+#ifdef CGAL_PMP_DEBUG_ORIENTATION
+      std::cout << "DEBUG: Set intersection is not empty\n";
+#endif
+      return false;
+    }
+
+    // set bit of compatible patches
+    for (F_cc_id id : bit_0_cc_set)
+    {
+      if (cc_handled[id])
+      {
+        if(cc_bits[id] == true)
+        {
+#ifdef CGAL_PMP_DEBUG_ORIENTATION
+          std::cout << "DEBUG: orientation bit already set to 1, incompatible with 0\n";
+#endif
+          return false;
+        }
+        else
+          continue;
+      }
+      cc_handled[id]=true;
+    }
+
+    // set bit of incompatible patches
+    for (F_cc_id id : bit_1_cc_set)
+    {
+      if (cc_handled[id])
+      {
+        if(cc_bits[id] == false)
+        {
+#ifdef CGAL_PMP_DEBUG_ORIENTATION
+          std::cout << "DEBUG: orientation bit already set to 0, incompatible with 1\n";
+#endif
+          return false;
+        }
+        else
+          continue;
+      }
+      cc_handled[id]=true;
+      cc_bits[id]=true;
+    }
+  }
+
+  // set the bit per face
+  for (face_descriptor f : faces(pm))
+    put(fbm, f, cc_bits[get(f_cc_ids,f)]);
+
+  return true;
+}
+
 } // namespace Polygon_mesh_processing
 } // namespace CGAL
 #endif // CGAL_ORIENT_POLYGON_MESH_H