diff --git a/Installation/CHANGES.md b/Installation/CHANGES.md
index 53170b00ebc..f720e98e1e3 100644
--- a/Installation/CHANGES.md
+++ b/Installation/CHANGES.md
@@ -1,5 +1,13 @@
 Release History
 ===============
+[Release 5.4](https://github.com/CGAL/cgal/releases/tag/v5.4)
+-----------
+
+Release date: December 2021
+
+### [Polygon Mesh Processing](https://doc.cgal.org/5.4/Manual/packages.html#PkgPolygonMeshProcessing)
+
+-   Added the function `CGAL::Polygon_mesh_processing::match_faces()`, which, given two polygon meshes, identifies their common faces as well as as faces present in only either of them.
 
 [Release 5.3](https://github.com/CGAL/cgal/releases/tag/v5.3)
 -----------
diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt
index 8e5e7c71448..9edc40b6e02 100644
--- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt
+++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt
@@ -198,6 +198,7 @@ The page \ref bgl_namedparameters "Named Parameters" describes their usage.
 - \link measure_grp `CGAL::Polygon_mesh_processing::edge_length()` \endlink
 - \link measure_grp `CGAL::Polygon_mesh_processing::face_border_length()` \endlink
 - \link measure_grp `CGAL::Polygon_mesh_processing::centroid()` \endlink
+- \link measure_grp `CGAL::Polygon_mesh_processing::match_faces()` \endlink
 
 \cgalCRPSection{Distance Functions}
 - `CGAL::Polygon_mesh_processing::approximate_Hausdorff_distance()`
diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt
index e853654b386..ad9c683900e 100644
--- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt
+++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt
@@ -90,6 +90,7 @@ create_single_source_cgal_program("locate_example.cpp")
 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")
 
 if(OpenMesh_FOUND)
 
diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/match_faces.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/match_faces.cpp
new file mode 100644
index 00000000000..9eeae6b6d00
--- /dev/null
+++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/match_faces.cpp
@@ -0,0 +1,58 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/Polyhedron_3.h>
+
+#include <CGAL/boost/graph/copy_face_graph.h>
+#include <CGAL/boost/graph/Euler_operations.h>
+#include <CGAL/Polygon_mesh_processing/measure.h>
+#include <CGAL/boost/graph/Named_function_parameters.h>
+#include <CGAL/boost/graph/named_params_helper.h>
+#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
+
+#include <fstream>
+#include <iostream>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel       K;
+
+typedef K::Point_3                                                Point;
+
+typedef CGAL::Surface_mesh<Point>                                 Surface_mesh;
+typedef CGAL::Polyhedron_3<K>                                     Polyhedron;
+typedef boost::graph_traits<Surface_mesh>::face_descriptor        face_descriptor_1;
+typedef boost::graph_traits<Polyhedron>::face_descriptor          face_descriptor_2;
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+int main(int argc, char* argv[])
+{
+  const char* filename1 = (argc > 1) ? argv[1] : "data/P.off";
+
+  Surface_mesh mesh1;
+  Polyhedron mesh2;
+  if(!PMP::IO::read_polygon_mesh(filename1, mesh1))
+  {
+    std::cerr << "Invalid input." << std::endl;
+    return 1;
+  }
+  CGAL::copy_face_graph(mesh1, mesh2);
+  CGAL::Euler::add_center_vertex(*halfedges(mesh2).begin(),mesh2);
+  std::vector<std::pair<face_descriptor_1, face_descriptor_2> > common;
+  std::vector<face_descriptor_1> m1_only;
+  std::vector<face_descriptor_2> m2_only;
+
+  PMP::match_faces(mesh1, mesh2, std::back_inserter(common), std::back_inserter(m1_only), std::back_inserter(m2_only));
+
+  std::cout <<"Faces only in m1 :"<< std::endl;
+  for(const face_descriptor_1& f : m1_only)
+    std::cout << " " << f;
+
+  std::cout <<"\n\nFaces only in m2:" << std::endl;
+  for(const face_descriptor_2& f : m2_only)
+    std::cout << " " << &(*f);
+
+  std::cout << "\n\nFaces in both:" << std::endl;
+  for(const auto& f_pair : common)
+    std::cout << " (" << f_pair.first << ", " << &(*f_pair.second);
+  std::cout << std::endl;
+
+  return 0;
+}
diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/measure.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/measure.h
index 911087ac7b9..f65caaaa87a 100644
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/measure.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/measure.h
@@ -28,10 +28,13 @@
 
 #include <CGAL/Lazy.h> // needed for CGAL::exact(FT)/CGAL::exact(Lazy_exact_nt<T>)
 
+#include <boost/container/small_vector.hpp>
 #include <boost/unordered_set.hpp>
 #include <boost/graph/graph_traits.hpp>
+#include <boost/dynamic_bitset.hpp>
 
 #include <utility>
+#include <algorithm>
 
 #ifdef DOXYGEN_RUNNING
 #define CGAL_PMP_NP_TEMPLATE_PARAMETERS NamedParameters
@@ -50,6 +53,14 @@ public:
 
 namespace Polygon_mesh_processing {
 
+namespace internal {
+
+inline void rearrange_face_ids(boost::container::small_vector<std::size_t, 4>& ids)
+{
+  auto min_elem = std::min_element(ids.begin(), ids.end());
+  std::rotate(ids.begin(), min_elem, ids.end());
+}
+}//namespace internal
 /**
   * \ingroup measure_grp
   * computes the length of an edge of a given polygon mesh.
@@ -820,6 +831,192 @@ centroid(const TriangleMesh& tmesh)
   return centroid(tmesh, CGAL::Polygon_mesh_processing::parameters::all_default());
 }
 
+
+/**
+  * \ingroup measure_grp
+  * identifies faces only present in `m1` and `m2` as well as the faces present
+  * in both polygon meshes. Two faces are matching if they have the same
+  * orientation and the same points.
+  *
+  * @tparam PolygonMesh1 a model of `HalfedgeListGraph` and `FaceListGraph`
+  * @tparam PolygonMesh2 a model of `HalfedgeListGraph` and `FaceListGraph`
+  * @tparam FaceOutputIterator1 model of `OutputIterator`
+     holding `boost::graph_traits<PolygonMesh1>::%face_descriptor`.
+  * @tparam FaceOutputIterator2 model of `OutputIterator`
+     holding `boost::graph_traits<PolygonMesh2>::%face_descriptor`.
+  * @tparam FacePairOutputIterator model of `OutputIterator`
+     holding `std::pair<boost::graph_traits<PolygonMesh1>::%face_descriptor,
+     boost::graph_traits<PolygonMesh2>::%face_descriptor`.
+  *
+  * @tparam NamedParameters1 a sequence of \ref bgl_namedparameters "Named Parameters"
+  * @tparam NamedParameters2 a sequence of \ref bgl_namedparameters "Named Parameters"
+  *
+  * @param m1 the first `PolygonMesh`
+  * @param m2 the second `PolygonMesh`
+  * @param common output iterator collecting the faces that are common to both meshes.
+  * @param m1_only output iterator collecting the faces that are only in `m1`
+  * @param m2_only output iterator collecting the faces that are only in `m2`
+  * @param np1 an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+  * @param np2 an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+  *
+  * \cgalNamedParamsBegin
+  *   \cgalParamNBegin{vertex_point_map}
+  *     \cgalParamDescription{a property map associating points to the vertices of `m1`}
+  *     \cgalParamType{a class model of `ReadablePropertyMap` with `boost::graph_traits<PolygonMesh1>::%vertex_descriptor`
+  *                    as key type and `%Point_3` as value type. `%Point_3` must be `LessThanComparable`.}
+  *     \cgalParamDefault{`boost::get(CGAL::vertex_point, m1)`}
+  *     \cgalParamExtra{The same holds for `m2` and `PolygonMesh2` and the point type must be the same for both meshes.}
+  *   \cgalParamNEnd
+  *
+  *   \cgalParamNBegin{vertex_index_map}
+  *     \cgalParamDescription{a property map associating to each vertex of `m1` a unique index between `0` and `num_vertices(m1) - 1`, and similarly for `m2`.}
+  *     \cgalParamType{a class model of `ReadablePropertyMap` with `boost::graph_traits<Graph>::%vertex_descriptor`
+  *                    as key type and `std::size_t` as value type}
+  *     \cgalParamDefault{an automatically indexed internal map}
+  *     \cgalParamExtra{If this parameter is not passed, internal machinery will create and initialize
+  *                     a face index property map, either using the internal property map if it exists
+  *                     or using an external map. The latter might result in  - slightly - worsened performance
+  *                     in case of non-constant complexity for index access. The same holds for `m2` and `PolygonMesh2`.}
+  *   \cgalParamNEnd
+  * \cgalNamedParamsEnd
+  *
+ */
+template< typename PolygonMesh1,
+          typename PolygonMesh2,
+          typename FacePairOutputIterator,
+          typename FaceOutputIterator1,
+          typename FaceOutputIterator2,
+          typename NamedParameters1,
+          typename NamedParameters2 >
+void match_faces(const PolygonMesh1& m1, const PolygonMesh2& m2,
+                 FacePairOutputIterator common, FaceOutputIterator1 m1_only, FaceOutputIterator2 m2_only,
+                 const NamedParameters1& np1, const NamedParameters2& np2)
+{
+  typedef typename GetVertexPointMap<PolygonMesh1, NamedParameters1>::const_type            VPMap1;
+  typedef typename GetVertexPointMap<PolygonMesh2, NamedParameters2>::const_type            VPMap2;
+  typedef typename GetInitializedVertexIndexMap<PolygonMesh1, NamedParameters1>::const_type VIMap1;
+  typedef typename GetInitializedVertexIndexMap<PolygonMesh2, NamedParameters2>::const_type VIMap2;
+  typedef typename boost::property_traits<VPMap2>::value_type                               Point_3;
+  typedef typename boost::graph_traits<PolygonMesh1>::face_descriptor                       face_descriptor_1;
+
+  using parameters::choose_parameter;
+  using parameters::get_parameter;
+
+  const VPMap1 vpm1 = choose_parameter(get_parameter(np1, internal_np::vertex_point),
+                                       get_const_property_map(vertex_point, m1));
+  const VPMap2 vpm2 = choose_parameter(get_parameter(np2, internal_np::vertex_point),
+                                       get_const_property_map(vertex_point, m2));
+  CGAL_static_assertion_msg((boost::is_same<typename boost::property_traits<VPMap1>::value_type,
+                             typename boost::property_traits<VPMap2>::value_type>::value),
+                            "Both vertex point maps must have the same point type.");
+
+  const VIMap1 vim1 = get_initialized_vertex_index_map(m1, np1);
+  const VIMap2 vim2 = get_initialized_vertex_index_map(m2, np2);
+
+  std::map<Point_3, std::size_t> point_id_map;
+
+  std::vector<std::size_t> m1_vertex_id(num_vertices(m1), -1);
+  std::vector<std::size_t> m2_vertex_id(num_vertices(m2), -1);
+  boost::dynamic_bitset<> shared_vertices(m1_vertex_id.size() + m2_vertex_id.size());
+
+  //iterate both meshes to set ids of all points, and set vertex/point_id maps.
+  std::size_t id = 0;
+  for(auto v : vertices(m1))
+  {
+    const typename boost::property_traits<VPMap1>::reference p = get(vpm1, v);
+    auto res = point_id_map.emplace(p, id);
+    if(res.second)
+      ++id;
+    m1_vertex_id[get(vim1, v)] = res.first->second;
+  }
+  for(auto v : vertices(m2))
+  {
+    const typename boost::property_traits<VPMap2>::reference p = get(vpm2, v);
+    auto res = point_id_map.emplace(p, id);
+    if(res.second)
+      ++id;
+    else
+      shared_vertices.set(res.first->second);
+    m2_vertex_id[get(vim2, v)] = res.first->second;
+  }
+
+  //fill a set with the "faces point-ids" of m1 and then iterate faces of m2 to compare.
+  std::map<boost::container::small_vector<std::size_t, 4>, face_descriptor_1> m1_faces_map;
+  for(auto f : faces(m1))
+  {
+    bool all_shared = true;
+    boost::container::small_vector<std::size_t, 4> ids;
+    for(auto v : CGAL::vertices_around_face(halfedge(f, m1), m1))
+    {
+      std::size_t vid = m1_vertex_id[get(vim1, v)];
+      ids.push_back(vid);
+      if(!shared_vertices.test(vid))
+      {
+        all_shared = false;
+        break;
+      }
+    }
+    if(all_shared)
+    {
+      internal::rearrange_face_ids(ids);
+      m1_faces_map.emplace(ids, f);
+    }
+    else
+      *m1_only++ = f;
+  }
+  for(auto f : faces(m2))
+  {
+    boost::container::small_vector<std::size_t, 4> ids;
+    bool all_shared = true;
+    for(auto v : CGAL::vertices_around_face(halfedge(f, m2), m2))
+    {
+      std::size_t vid = m2_vertex_id[get(vim2, v)];
+      ids.push_back(vid);
+      if(!shared_vertices.test(vid))
+      {
+        all_shared = false;
+        break;
+      }
+    }
+    if(all_shared)
+    {
+      internal::rearrange_face_ids(ids);
+      auto it = m1_faces_map.find(ids);
+      if(it != m1_faces_map.end())
+      {
+        *common++ = std::make_pair(it->second, f);
+        m1_faces_map.erase(it);
+      }
+      else
+      {
+        *m2_only++ = f;
+      }
+    }
+    else
+      *m2_only++ = f;
+  }
+  //all shared faces have been removed from the map, so all that remains must go in m1_only
+  for(const auto& it : m1_faces_map)
+  {
+    *m1_only++ = it.second;
+  }
+}
+
+template<typename PolygonMesh1, typename PolygonMesh2, typename FacePairOutputIterator, typename FaceOutputIterator1, typename FaceOutputIterator2, typename NamedParameters>
+void match_faces(const PolygonMesh1& m1, const PolygonMesh2& m2,
+                 FacePairOutputIterator common, FaceOutputIterator1 m1_only, FaceOutputIterator2 m2_only,
+                 const NamedParameters& np)
+{
+  match_faces(m1, m2, common, m1_only, m2_only, np, parameters::all_default());
+}
+
+template<typename PolygonMesh1, typename PolygonMesh2, typename FacePairOutputIterator, typename FaceOutputIterator1, typename FaceOutputIterator2>
+void match_faces(const PolygonMesh1& m1, const PolygonMesh2& m2,
+                 FacePairOutputIterator common, FaceOutputIterator1 m1_only, FaceOutputIterator2 m2_only)
+{
+  match_faces(m1, m2, common, m1_only, m2_only, parameters::all_default(), parameters::all_default());
+}
+
 } // namespace Polygon_mesh_processing
 } // namespace CGAL
 
diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/data/cube_quad2.off b/Polygon_mesh_processing/test/Polygon_mesh_processing/data/cube_quad2.off
new file mode 100644
index 00000000000..45cb8852a90
--- /dev/null
+++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/data/cube_quad2.off
@@ -0,0 +1,18 @@
+OFF
+8 6 0
+
+-1 -1 -1
+-1 1 -1
+1 1 -1
+1 -1 -1
+-1.53485 -0.408879 0.387354
+-1 1 1
+1 1 1
+1 -1 1
+
+4  0 3 7 4
+4  3 2 6 7
+4  2 1 5 6
+4  1 0 4 5
+4  4 7 6 5
+4  0 1 2 3
diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/data/tri1-hole.off b/Polygon_mesh_processing/test/Polygon_mesh_processing/data/tri1-hole.off
new file mode 100644
index 00000000000..7ba3bb6021b
--- /dev/null
+++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/data/tri1-hole.off
@@ -0,0 +1,37 @@
+OFF
+16 17 0
+
+0 0 0
+0.333333 0 0
+0.666667 0 0
+1 0 0
+0 0.333333 0.333333
+0.333333 0.333333 0.333333
+0.666667 0.333333 0.333333
+1 0.333333 0.333333
+0 0.666667 0.666667
+0.333333 0.666667 0.666667
+0.666667 0.666667 0.666667
+1 0.666667 0.666667
+0 1 1
+0.333333 1 1
+0.666667 1 1
+1 1 1
+
+3  0 1 4
+3  1 5 4
+3  4 5 8
+3  5 9 8
+3  8 9 12
+3  9 13 12
+3  2 6 5
+3  5 6 9
+3  6 10 9
+3  9 10 13
+3  10 14 13
+3  2 3 6
+3  3 7 6
+3  6 7 10
+3  7 11 10
+3  10 11 14
+3  11 15 14
diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/data/tri1.off b/Polygon_mesh_processing/test/Polygon_mesh_processing/data/tri1.off
new file mode 100644
index 00000000000..2aada48783f
--- /dev/null
+++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/data/tri1.off
@@ -0,0 +1,38 @@
+OFF
+16 18 0
+
+0 0 0
+0.333333 0 0
+0.666667 0 0
+1 0 0
+0 0.333333 0.333333
+0.333333 0.333333 0.333333
+0.666667 0.333333 0.333333
+1 0.333333 0.333333
+0 0.666667 0.666667
+0.333333 0.666667 0.666667
+0.666667 0.666667 0.666667
+1 0.666667 0.666667
+0 1 1
+0.333333 1 1
+0.666667 1 1
+1 1 1
+
+3  0 1 4
+3  1 5 4
+3  4 5 8
+3  5 9 8
+3  8 9 12
+3  9 13 12
+3  1 2 5
+3  2 6 5
+3  5 6 9
+3  6 10 9
+3  9 10 13
+3  10 14 13
+3  2 3 6
+3  3 7 6
+3  6 7 10
+3  7 11 10
+3  10 11 14
+3  11 15 14
diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/data/tri2-out.off b/Polygon_mesh_processing/test/Polygon_mesh_processing/data/tri2-out.off
new file mode 100644
index 00000000000..1127abddf90
--- /dev/null
+++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/data/tri2-out.off
@@ -0,0 +1,38 @@
+OFF
+16 18 0
+
+0 0 0
+0.33333299999999999 0 0
+0.66666700000000001 0 0
+1 0 0
+0 0.33333299999999999 0.33333299999999999
+0.33333299999999999 0.33333299999999999 0.33333299999999999
+0.67039400000000005 0.27928500000000001 0.33333299999999999
+1 0.33333299999999999 0.33333299999999999
+0 0.66666700000000001 0.66666700000000001
+0.33333299999999999 0.66666700000000001 0.66666700000000001
+0.66666700000000001 0.66666700000000001 0.66666700000000001
+1 0.66666700000000001 0.66666700000000001
+0 1 1
+0.37796099999999999 1.08385 1
+0.73563699999999999 1.01488 0.94896599999999998
+0.96348599999999995 1.0906100000000001 1
+
+3  4 1 0
+3  4 5 1
+3  8 5 4
+3  8 9 5
+3  12 9 8
+3  12 13 9
+3  5 2 1
+3  5 6 2
+3  9 6 5
+3  9 10 6
+3  13 10 9
+3  13 14 10
+3  6 3 2
+3  6 7 3
+3  10 7 6
+3  10 11 7
+3  14 11 10
+3  14 15 11
diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/data/tri2.off b/Polygon_mesh_processing/test/Polygon_mesh_processing/data/tri2.off
new file mode 100644
index 00000000000..462ea37293a
--- /dev/null
+++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/data/tri2.off
@@ -0,0 +1,38 @@
+OFF
+16 18 0
+
+0 0 0
+0.333333 0 0
+0.666667 0 0
+1 0 0
+0 0.333333 0.333333
+0.333333 0.333333 0.333333
+0.670394 0.279285 0.333333
+1 0.333333 0.333333
+0 0.666667 0.666667
+0.333333 0.666667 0.666667
+0.666667 0.666667 0.666667
+1 0.666667 0.666667
+0 1 1
+0.377961 1.08385 1
+0.735637 1.01488 0.948966
+0.963486 1.09061 1
+
+3  0 1 4
+3  1 5 4
+3  4 5 8
+3  5 9 8
+3  8 9 12
+3  9 13 12
+3  1 2 5
+3  2 6 5
+3  5 6 9
+3  6 10 9
+3  9 10 13
+3  10 14 13
+3  2 3 6
+3  3 7 6
+3  6 7 10
+3  7 11 10
+3  10 11 14
+3  11 15 14
diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/measures_test.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/measures_test.cpp
index b1537349cc1..b8db4721474 100644
--- a/Polygon_mesh_processing/test/Polygon_mesh_processing/measures_test.cpp
+++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/measures_test.cpp
@@ -203,6 +203,150 @@ void test_centroid(const char* filename)
 
 }
 
+template <typename PolygonMesh1, typename PolygonMesh2 >
+void test_compare()
+{
+  typedef typename boost::graph_traits<PolygonMesh1>::face_descriptor        face_descriptor1;
+  typedef typename boost::graph_traits<PolygonMesh2>::face_descriptor        face_descriptor2;
+  namespace PMP = CGAL::Polygon_mesh_processing;
+
+  PolygonMesh1 mesh1;
+  PolygonMesh2 mesh2;
+  std::vector<std::pair<face_descriptor1, face_descriptor2> > common;
+  common.clear();
+  std::vector<face_descriptor1> m1_only;
+  std::vector<face_descriptor2> m2_only;
+  /*************************
+   * triangulated and open *
+   * **********************/
+
+  std::ifstream input("data/tri1.off");
+  if(! (input >> mesh1))
+  {
+    std::cerr << "Invalid input." << std::endl;
+    assert (false);
+    return;
+  }
+  input.close();
+  input.open("data/tri2.off");
+  if(! (input >> mesh2))
+  {
+    std::cerr << "Invalid input." << std::endl;
+    assert (false);
+    return;
+  }
+  input.close();
+  PMP::match_faces(mesh1, mesh2, std::back_inserter(common), std::back_inserter(m1_only), std::back_inserter(m2_only), CGAL::parameters::all_default(), CGAL::parameters::all_default());
+  assert(common.size() == 7);
+  assert(m1_only.size() == 11);
+  assert(m2_only.size() == 11);
+  /*************************
+   **** quad and closed ****
+   * **********************/
+  CGAL::clear(mesh1);
+  CGAL::clear(mesh2);
+  common.clear();
+  m1_only.clear();
+  m2_only.clear();
+  input.open("data/cube_quad.off");
+  if(! (input >> mesh1))
+  {
+    std::cerr << "Invalid input." << std::endl;
+    assert (false);
+    return;
+  }
+  input.close();
+  input.open("data/cube_quad2.off");
+  if(! (input >> mesh2))
+  {
+    std::cerr << "Invalid input." << std::endl;
+    assert (false);
+    return;
+  }
+  input.close();
+  std::unordered_map<face_descriptor1, std::size_t> fim1;
+  std::unordered_map<face_descriptor2, std::size_t> fim2;
+  std::size_t id = 0;
+  for(const auto& f : faces(mesh1))
+  {
+    fim1.insert(std::make_pair(f, id++));
+  }
+  id = 0;
+  for(const auto& f : faces(mesh2))
+  {
+    fim2.insert(std::make_pair(f, id++));
+  }
+  PMP::match_faces(mesh1, mesh2, std::back_inserter(common), std::back_inserter(m1_only), std::back_inserter(m2_only), CGAL::parameters::all_default(), CGAL::parameters::all_default());
+  assert(common.size() == 3);
+  assert(m1_only.size() == 3);
+  assert(fim1[m1_only[0]] == 0);
+  assert(fim1[m1_only[1]] == 3);
+  assert(fim1[m1_only[2]] == 4);
+  assert(m2_only.size() == 3);
+  assert(fim2[m2_only[0]] == 0);
+  assert(fim2[m2_only[1]] == 3);
+  assert(fim2[m2_only[2]] == 4);
+  /*************************
+   **** tri and hole****
+   * **********************/
+  CGAL::clear(mesh1);
+  CGAL::clear(mesh2);
+  common.clear();
+  m1_only.clear();
+  m2_only.clear();
+  input.open("data/tri1.off");
+  if(! (input >> mesh1))
+  {
+    std::cerr << "Invalid input." << std::endl;
+    assert (false);
+    return;
+  }
+  input.close();
+  input.open("data/tri1-hole.off");
+  if(! (input >> mesh2))
+  {
+    std::cerr << "Invalid input." << std::endl;
+    assert (false);
+    return;
+  }
+  input.close();
+  PMP::match_faces(mesh1, mesh2, std::back_inserter(common), std::back_inserter(m1_only), std::back_inserter(m2_only), CGAL::parameters::all_default(), CGAL::parameters::all_default());
+  assert(common.size() == 17);
+  assert(m1_only.size() == 1);
+  assert(m2_only.size() == 0);
+
+  /*************************
+   **** tri and orient****
+   * **********************/
+  CGAL::clear(mesh1);
+  CGAL::clear(mesh2);
+  common.clear();
+  m1_only.clear();
+  m2_only.clear();
+  input.open("data/tri2.off");
+  if(! (input >> mesh1))
+  {
+    std::cerr << "Invalid input." << std::endl;
+    assert (false);
+    return;
+  }
+  input.close();
+  input.open("data/tri2-out.off");
+  if(! (input >> mesh2))
+  {
+    std::cerr << "Invalid input." << std::endl;
+    assert (false);
+    return;
+  }
+  input.close();
+  PMP::match_faces(mesh1, mesh2, std::back_inserter(common), std::back_inserter(m1_only),
+                      std::back_inserter(m2_only));
+  assert(common.size() == 0);
+  assert(m1_only.size() == 18);
+  assert(m2_only.size() == 18);
+
+}
+
 int main(int argc, char* argv[])
 {
   const char* filename_polyhedron =
@@ -218,7 +362,11 @@ int main(int argc, char* argv[])
   // It won't work with Epec for large meshes as it builds up a deep DAG
   // leading to a stackoverflow when the destructor is called.
   test_centroid<CGAL::Surface_mesh<Epic::Point_3>,Epic>(filename_surface_mesh);
-
+  test_compare<CGAL::Polyhedron_3<Epic>, CGAL::Surface_mesh<Epic::Point_3> >();
+  test_compare<CGAL::Polyhedron_3<Epec>, CGAL::Surface_mesh<Epec::Point_3> >();
+  test_compare<CGAL::Surface_mesh<Epic::Point_3>, CGAL::Polyhedron_3<Epic> >();
+  test_compare<CGAL::Surface_mesh<Epec::Point_3>, CGAL::Polyhedron_3<Epec> >();
   std::cerr << "All done." << std::endl;
   return 0;
 }
+
diff --git a/Polyhedron/demo/Polyhedron/Plugins/Operations_on_polyhedra/Diff_between_meshes_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/Operations_on_polyhedra/Diff_between_meshes_plugin.cpp
index 3147fca83cf..8ab19827c47 100644
--- a/Polyhedron/demo/Polyhedron/Plugins/Operations_on_polyhedra/Diff_between_meshes_plugin.cpp
+++ b/Polyhedron/demo/Polyhedron/Plugins/Operations_on_polyhedra/Diff_between_meshes_plugin.cpp
@@ -8,10 +8,12 @@
 #include "Scene_surface_mesh_item.h"
 
 #include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
+#include <CGAL/Polygon_mesh_processing/measure.h>
 
 #include <CGAL/Three/Polyhedron_demo_plugin_interface.h>
 #include <CGAL/Three/Scene_group_item.h>
 #include <CGAL/Three/Three.h>
+#include <CGAL/boost/graph/Face_filtered_graph.h>
 #include "Messages_interface.h"
 
 using namespace CGAL::Three;
@@ -48,7 +50,6 @@ public Q_SLOTS:
 private:
   CGAL::Three::Scene_interface* scene;
   QAction* actionDiff;
-  SMesh* diff(SMesh* m1, SMesh* m2, bool compute_common);
 
 }; // end Polyhedron_demo_diff_between_meshes_plugin
 
@@ -67,95 +68,12 @@ QList<QAction*> Polyhedron_demo_diff_between_meshes_plugin::actions() const {
   return QList<QAction*>() << actionDiff;
 }
 
-SMesh* Polyhedron_demo_diff_between_meshes_plugin::diff(SMesh* m1, SMesh* m2, bool compute_common = false)
-{
-  std::map<Point_3, std::size_t> point_id_map;
-  std::vector<std::size_t> m1_vertex_id(num_vertices(*m1), -1);
-  std::vector<std::size_t> m2_vertex_id(num_vertices(*m2), -1);
-
-  //iterate both meshes to set ids to all points, and set vertex/point_id maps.
-  std::size_t id =0;
-  for(auto v : m1->vertices())
-  {
-    Point_3 p = m1->point(v);
-    auto res = point_id_map.insert(std::make_pair(p, id));
-    if(res.second)
-      id++;
-    m1_vertex_id[(std::size_t)v]=res.first->second;
-  }
-  for(auto v : m2->vertices())
-  {
-    Point_3 p = m2->point(v);
-    auto res = point_id_map.insert(std::make_pair(p, id));
-    if(res.second)
-      id++;
-    m2_vertex_id[(std::size_t)v]=res.first->second;
-  }
-
-  //fill a set with the "faces point-ids" of m1 and then iterate faces of m2 to compare.
-  std::set<std::vector<std::size_t> > m1_faces;
-  for(auto f : m1->faces())
-  {
-    std::vector<std::size_t> ids;
-    for(auto v : CGAL::vertices_around_face(halfedge(f, *m1), *m1))
-    {
-      ids.push_back(m1_vertex_id[(std::size_t)v]);
-    }
-    std::sort(ids.begin(), ids.end());
-    m1_faces.insert(ids);
-  }
-
-  std::vector<SMesh::Face_index> common_faces;
-  std::vector<Point_3> common_points;
-  std::map<SMesh::Vertex_index, std::size_t> id_map;
-  id = 0;
-
-  for(auto f : m2->faces())
-  {
-    std::vector<std::size_t> ids;
-    for(auto v : CGAL::vertices_around_face(halfedge(f, *m2), *m2))
-    {
-      ids.push_back(m2_vertex_id[(std::size_t)v]);
-    }
-    std::sort(ids.begin(), ids.end());
-    if(!((m1_faces.find(ids) != m1_faces.end()) ^ compute_common))
-    {
-      common_faces.push_back(f);
-      for(auto v : CGAL::vertices_around_face(halfedge(f, *m2), *m2))
-      {
-        auto res = id_map.insert(std::make_pair(v,id));
-        if(res.second)
-        {
-          common_points.push_back(m2->point(v));
-          id++;
-        }
-      }
-    }
-  }
-
-  //iterate m1_faces and fill a polygon vector using the id_map previously filled.
-  std::vector<std::vector<std::size_t> > polygons(common_faces.size());
-  id = 0;
-  for(auto f : common_faces)
-  {
-    for(auto v : vertices_around_face(halfedge(f, *m2),*m2))
-    {
-      polygons[id].push_back(id_map[v]);
-    }
-    ++id;
-  }
-
-  SMesh* common = new SMesh();
-  CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh<SMesh>(
-    common_points, polygons, *common);
-  return common;
-
-}
-
 
 void Polyhedron_demo_diff_between_meshes_plugin::diff()
 {
 
+  typedef CGAL::Face_filtered_graph<SMesh> Filtered_graph;
+
   QCursor c(Qt::WaitCursor);
   CGAL::Three::Three::CursorScopeGuard guard(c);
 
@@ -166,32 +84,64 @@ void Polyhedron_demo_diff_between_meshes_plugin::diff()
       *m2_item = qobject_cast<Scene_surface_mesh_item*>(
         scene->item(scene->selectionIndices().back()));
 
-  SMesh* m1=m1_item->face_graph(),
-      *m2=m2_item->face_graph();
+  SMesh m1=*m1_item->face_graph(),
+      m2=*m2_item->face_graph();
+  std::vector<face_descriptor> m1_only, m2_only;
+  std::vector<std::pair<face_descriptor, face_descriptor> > common;
+  CGAL::Polygon_mesh_processing::match_faces(
+        m1,
+        m2,
+        std::back_inserter(common),
+        std::back_inserter(m1_only),
+        std::back_inserter(m2_only));
 
-  SMesh* m1_over_m2 = diff(m1, m2);
-  SMesh* m2_over_m1 = diff(m2, m1);
-  SMesh* common = diff(m2, m1, true);
-
-  Scene_surface_mesh_item* m1_over_m2_item = new Scene_surface_mesh_item(m1_over_m2);
-  m1_over_m2_item->setColor(QColor(Qt::blue));
-  m1_over_m2_item->setName(QString("%2 - %1").arg(m1_item->name()).arg(m2_item->name()));
-  CGAL::Three::Three::scene()->addItem(m1_over_m2_item);
-  Scene_surface_mesh_item* m2_over_m1_item = new Scene_surface_mesh_item(m2_over_m1);
-  m2_over_m1_item->setColor(QColor(Qt::red));
-  m2_over_m1_item->setName(QString("%1 - %2").arg(m1_item->name()).arg(m2_item->name()));
-  CGAL::Three::Three::scene()->addItem(m2_over_m1_item);
-  Scene_surface_mesh_item* common_item = new Scene_surface_mesh_item(common);
-  common_item->setColor(QColor(Qt::green));
-  CGAL::Three::Three::scene()->addItem(common_item);
-  common_item->setName(QString("%1 && %2").arg(m1_item->name()).arg(m2_item->name()));
+  Filtered_graph filter1(m1, m1_only);
+  SMesh mesh1_only, mesh2_only, common_mesh;
+  CGAL::copy_face_graph(filter1, mesh1_only);
+  Scene_surface_mesh_item* mesh1_only_item = nullptr;
+  if(mesh1_only.faces().size() > 0)
+  {
+    mesh1_only_item = new Scene_surface_mesh_item(mesh1_only);
+    mesh1_only_item->setColor(QColor(Qt::blue));
+    mesh1_only_item->setName(QString("%1_only").arg(m1_item->name()));
+    CGAL::Three::Three::scene()->addItem(mesh1_only_item);
+  }
 
+  Filtered_graph filter2(m2, m2_only);
+  CGAL::copy_face_graph(filter2, mesh2_only);
+  Scene_surface_mesh_item* mesh2_only_item = nullptr;
+  if(mesh2_only.faces().size() > 0)
+  {
+    mesh2_only_item = new Scene_surface_mesh_item(mesh2_only);
+    mesh2_only_item->setColor(QColor(Qt::red));
+    mesh2_only_item->setName(QString("%1_only").arg(m2_item->name()));
+    CGAL::Three::Three::scene()->addItem(mesh2_only_item);
+  }
+  m1_only.clear();
+  m1_only.reserve(common.size());
+  for(const auto& f_pair : common)
+  {
+    m1_only.push_back(f_pair.first);
+  }
+  Filtered_graph filter_common(m1, m1_only);
+  CGAL::copy_face_graph(filter_common, common_mesh);
+  Scene_surface_mesh_item* common_item = nullptr;
+  if(common_mesh.faces().size() > 0)
+  {
+    common_item = new Scene_surface_mesh_item(common_mesh);
+    common_item->setColor(QColor(Qt::green));
+    CGAL::Three::Three::scene()->addItem(common_item);
+    common_item->setName(QString("%1 && %2").arg(m1_item->name()).arg(m2_item->name()));
+  }
   Scene_group_item* group = new Scene_group_item();
   group->setName("Diff result");
   CGAL::Three::Three::scene()->addItem(group);
-  CGAL::Three::Three::scene()->changeGroup(m1_over_m2_item, group);
-  CGAL::Three::Three::scene()->changeGroup(m2_over_m1_item, group);
-  CGAL::Three::Three::scene()->changeGroup(common_item, group);
+  if(mesh1_only_item)
+    CGAL::Three::Three::scene()->changeGroup(mesh1_only_item, group);
+  if(mesh2_only_item)
+    CGAL::Three::Three::scene()->changeGroup(mesh2_only_item, group);
+  if(common_item)
+    CGAL::Three::Three::scene()->changeGroup(common_item, group);
 
   m1_item->setVisible(false);
   m2_item->setVisible(false);