PMP: compare_meshes

Installation/CHANGES.md

@@ -1,5 +1,14 @@
 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()` that given two polygon meshes identifies faces present in only one of the two meshes as well as faces present in both.
+
 
 [Release 5.3](https://github.com/CGAL/cgal/releases/tag/v5.3)
 -----------
@@ -459,7 +468,7 @@ Release date: November 2019
   - LLVM Clang version 8.0 or later (on Linux or MacOS), and
   - Apple Clang compiler versions 7.0.2 and 10.0.1 (on MacOS).
 - Since CGAL 4.9, CGAL can be used as a header-only library, with
-  dependencies. Since CGAL 5.0, that is now the default, unless
+  dependencies. Since CGAL 5.0, that is now the default, unless
   specified differently in the (optional) CMake configuration.
 - The section "Getting Started with CGAL" of the documentation has
   been updated and reorganized.
@@ -2299,14 +2308,14 @@ Release date: October 2014
 -   Changes in the set of supported platforms:
     -   The Microsoft Windows Visual C++ compiler 2008 (VC9) is no
         longer supported since CGAL-4.5.
--   Since CGAL version 4.0, Eigen was the recommended third-party
+-   Since CGAL version 4.0, Eigen was the recommended third-party
     library to use with *Planar Parameterization of Triangulated Surface
     Meshes*, *Surface Reconstruction from Point Sets*, *Approximation of
     Ridges and Umbilics on Triangulated Surface Meshes*, and *Estimation
     of Local Differential Properties of Point-Sampled Surfaces*
-    packages. From CGAL version 4.5, Taucs, Blas and Lapack are no
+    packages. From CGAL version 4.5, Taucs, Blas and Lapack are no
     longer supported.
--   CGAL is now compatible with the new CMake version 3.0.
+-   CGAL is now compatible with the new CMake version 3.0.
 
 ### Triangulated Surface Mesh Deformation (new package)
 
@@ -2445,7 +2454,7 @@ Release date: April 2014
 
 -   Additional supported platforms:
     -   The Apple Clang compiler version 5.0 is now supported on
-        OS X Mavericks.
+        OS X Mavericks.
     -   The Microsoft Windows Visual C++ compiler 2013 (VC12) is now
         supported.
 
@@ -2583,7 +2592,7 @@ Release date: October 2013
     transparent to the user thanks to the implicit constructor added to
     `CGAL::Object`. However, it is recommended to upgrade your code. The
     previous behavior can be restored by defining the macro
-    `CGAL_INTERSECTION_VERSION` to 1.
+    `CGAL_INTERSECTION_VERSION` to 1.
 
 #### 2D Arrangements
 
@@ -2844,7 +2853,7 @@ Release date: October 2012
 
 -   Additional supported platforms:
     -   The Apple Clang compiler versions 3.1 and 3.2 are now supported
-        on Mac OS X.
+        on Mac OS X.
 -   Improved configuration for essential and optional external third
     party software
 -   Added more general script to create CMakeLists.txt files:

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()`

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("compare_meshes_example.cpp")
 
 if(OpenMesh_FOUND)
 

Polygon_mesh_processing/examples/Polygon_mesh_processing/compare_meshes_example.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 auto& f : m1_only)
+    std::cout << " " << f;
+
+  std::cout<<"\n\nFaces only in m2:"<<std::endl;
+  for(const auto& 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;
+}

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 @@ class GetGeomTraits<CGAL_PMP_NP_CLASS, NP>
 
 namespace Polygon_mesh_processing {
 
+namespace pmp_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());
+}
+}//end pmp_internal
 /**
   * \ingroup measure_grp
   * computes the length of an edge of a given polygon mesh.
@@ -820,6 +831,186 @@ 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 identical 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 `PolygonMesh1`
+  * @param m2 the second `PolygonMesh2`
+  * @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 FaceOutputIterator1,
+          typename FaceOutputIterator2,
+          typename FacePairOutputIterator,
+          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)
+{
+  using parameters::choose_parameter;
+  using parameters::get_parameter;
+  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;
+  VPMap1 vpm1 = choose_parameter(get_parameter(np1, internal_np::vertex_point),
+                                      get_const_property_map(vertex_point, m1));
+  VPMap2 vpm2 = choose_parameter(get_parameter(np2, internal_np::vertex_point),
+                                      get_const_property_map(vertex_point, m2));
+  VIMap1 vim1 = get_initialized_vertex_index_map(m1, np1);
+  VIMap2 vim2 = get_initialized_vertex_index_map(m2, np2);
+  typedef typename boost::property_traits<VPMap2>::value_type Point_3;
+  typedef typename boost::graph_traits<PolygonMesh1>::face_descriptor face_descriptor_1;
+
+  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 Point_3& p = get(vpm1, v);
+    auto res = point_id_map.insert(std::make_pair(p, id));
+    if(res.second)
+      ++id;
+    m1_vertex_id[get(vim1, v)]=res.first->second;
+  }
+  for(auto v : vertices(m2))
+  {
+    const Point_3& p = get(vpm2, v);
+    auto res = point_id_map.insert(std::make_pair(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)
+    {
+      pmp_internal::rearrange_face_ids(ids);
+      m1_faces_map.insert({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)
+    {
+      pmp_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 FaceOutputIterator1, typename FaceOutputIterator2, typename FacePairOutputIterator, 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 FaceOutputIterator1, typename FaceOutputIterator2, typename FacePairOutputIterator>
+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
 

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