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::compare_meshes()` which allows, given two meshes, to geometrically detect which faces are only in the first, only in the second, or in both meshes.
+
 
 [Release 5.3](https://github.com/CGAL/cgal/releases/tag/v5.3)
 -----------

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,56 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Surface_mesh.h>
+
+#include <CGAL/Polygon_mesh_processing/measure.h>
+#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
+#include <CGAL/boost/graph/Named_function_parameters.h>
+#include <CGAL/boost/graph/named_params_helper.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 boost::graph_traits<Surface_mesh>::vertex_descriptor      vertex_descriptor;
+typedef boost::graph_traits<Surface_mesh>::face_descriptor        face_descriptor;
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+int main(int argc, char* argv[])
+{
+  const char* filename1 = (argc > 1) ? argv[1] : "data/tet.off";
+  const char* filename2 = (argc > 2) ? argv[2] : "data/tet-bis.off";
+
+  Surface_mesh mesh1, mesh2;
+  if(!PMP::read_polygon_mesh(filename1, mesh1))
+  {
+    std::cerr << "Invalid input." << std::endl;
+    return 1;
+  }
+  if(!PMP::read_polygon_mesh(filename2, mesh2))
+  {
+    std::cerr << "Invalid input." << std::endl;
+    return 1;
+  }
+  std::vector<std::pair<face_descriptor, face_descriptor> > common;
+  std::vector<face_descriptor> m1_only, m2_only;
+  PMP::compare_meshes(mesh1, mesh2, std::back_inserter(common), std::back_inserter(m1_only), std::back_inserter(m2_only), CGAL::parameters::all_default(), CGAL::parameters::all_default());
+  std::cout<<"Faces only in m1 : "<<std::endl;
+  for(const auto& f : m1_only)
+  {
+    std::cout<<f<<", ";
+  }
+  std::cout<<"\n Faces only in m2: "<<std::endl;
+  for(const auto& f : m2_only)
+  {
+    std::cout<<f<<", ";
+  }
+  std::cout<<"\n Faces in both: "<<std::endl;
+  for(const auto& f_pair : common)
+  {
+    std::cout<<f_pair.first<<", "<<f_pair.second<<";;"<<std::endl;
+  }
+  return 0;
+}

Polygon_mesh_processing/examples/Polygon_mesh_processing/data/tet-bis.off

@@ -0,0 +1,12 @@
+OFF
+4 4 0
+
+0 0 0
+1 0 0
+0 1 0
+0 0 1
+
+3  2 1 0
+3  0 3 2
+3  2 3 1
+3  1 3 0

Polygon_mesh_processing/examples/Polygon_mesh_processing/data/tet.off

@@ -0,0 +1,12 @@
+OFF
+4 4 0
+
+0.0521823 0.0693697 0.0558329
+1 0 0
+0 1 0
+0 0 1
+
+3  2 1 0
+3  0 3 2
+3  2 3 1
+3  1 3 0

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/measure.h

@@ -28,8 +28,10 @@
 
 #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>
 
@@ -820,6 +822,156 @@ centroid(const TriangleMesh& tmesh)
   return centroid(tmesh, CGAL::Polygon_mesh_processing::parameters::all_default());
 }
 
+
+/**
+  * \ingroup measure_grp
+  * given two meshes, separates the faces that are only in one, the faces
+  * that are only in the other one, and the faces that are common to both.
+  * The orientation of the faces is ignored during the comparison.
+  *
+  * @tparam PolygonMesh a model of `HalfedgeListGraph` and `FaceListGraph`
+  * @tparam OutputFaceIterator model of `OutputIterator`
+     holding `boost::graph_traits<PolygonMesh>::%face_descriptor`
+     for faces that are only in one mesh.
+  * @tparam OutputFacePairIterator model of `OutputIterator`
+     holding `std::pair<boost::graph_traits<PolygonMesh>::%face_descriptor,
+     boost::graph_traits<PolygonMesh>::%face_descriptor`
+     for faces that are shared by both meshes.
+  *
+  * @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` (`m2`)}
+  *     \cgalParamType{a class model of `ReadablePropertyMap` with `boost::graph_traits<PolygonMesh>::%vertex_descriptor`
+  *                    as key type and `%Point_3` as value type. `%Point_3` must be LessThanComparable.}
+  *     \cgalParamDefault{`boost::get(CGAL::vertex_point, m1 (m2))`}
+  *   \cgalParamNEnd
+  *
+  *   \cgalParamNBegin{vertex_index_map}
+  *     \cgalParamDescription{a property map associating to each vertex of `m1` (`m2`) a unique index between `0` and `num_vertices(m1 (m2)) - 1`}
+  *     \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.}
+  *   \cgalParamNEnd
+  * \cgalNamedParamsEnd
+  *
+ */
+template<typename PolygonMesh, typename OutputFaceIterator, typename OutputFacePairIterator, typename NamedParameters1, typename NamedParameters2 >
+void compare_meshes(const PolygonMesh& m1, const PolygonMesh& m2,
+                    OutputFacePairIterator common, OutputFaceIterator m1_only, OutputFaceIterator m2_only,
+                    const NamedParameters1& np1,const NamedParameters2& np2)
+{
+  using parameters::choose_parameter;
+  using parameters::get_parameter;
+  typedef typename GetVertexPointMap < PolygonMesh, NamedParameters1>::const_type VPMap1;
+  typedef typename GetVertexPointMap < PolygonMesh, NamedParameters2>::const_type VPMap2;
+  typedef typename GetInitializedVertexIndexMap<PolygonMesh, NamedParameters1>::const_type VIMap1;
+  typedef typename GetInitializedVertexIndexMap<PolygonMesh, 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);
+  VIMap1 vim2 = get_initialized_vertex_index_map(m2, np2);
+  typedef typename boost::property_traits<VPMap2>::value_type Point_3;
+  typedef typename boost::graph_traits<PolygonMesh>::face_descriptor face_descriptor;
+
+  std::map<Point_3, std::size_t> point_id_map;
+
+  boost::dynamic_bitset<> shared_vertices(num_vertices(m1) + num_vertices(m2));
+  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 : 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[static_cast<std::size_t>(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[static_cast<std::size_t>(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> 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[static_cast<std::size_t>(get(vim1, v))];
+      ids.push_back(vid);
+      if(!shared_vertices.test(vid))
+        all_shared = false;
+    }
+    std::sort(ids.begin(), ids.end());
+    if(all_shared)
+      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[static_cast<std::size_t>(get(vim2, v))];
+      ids.push_back(vid);
+      if(!shared_vertices.test(vid))
+      {
+        all_shared = false;
+      }
+    }
+    std::sort(ids.begin(), ids.end());
+    if(all_shared)
+      *common++ = std::make_pair(m1_faces_map[ids], f);
+    else
+      *m2_only++ = f;
+  }
+}
+
+template<typename PolygonMesh, typename OutputFaceIterator, typename OutputFacePairIterator, typename NamedParameters>
+void compare_meshes(const PolygonMesh& m1, const PolygonMesh& m2,
+                    OutputFacePairIterator common, OutputFaceIterator m1_only, OutputFaceIterator m2_only,
+                    const NamedParameters& np)
+{
+  compare_meshes(m1, m2, common, m1_only, m2_only, np, parameters::all_default());
+}
+
+template<typename PolygonMesh, typename OutputFaceIterator, typename OutputFacePairIterator>
+void compare_meshes(const PolygonMesh& m1, const PolygonMesh& m2,
+                    OutputFacePairIterator common, OutputFaceIterator m1_only, OutputFaceIterator m2_only)
+{
+  compare_meshes(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

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

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