first version of keep_largest_connected_component

BGL/include/CGAL/boost/graph/Dual.h

@@ -78,6 +78,13 @@ struct property_map<CGAL::Dual<P>, boost::vertex_index_t>
   typedef typename property_map<P, CGAL::face_index_t>::const_type const_type; 
 };
 
+template <typename P>
+struct property_map<CGAL::Dual<P>, boost::face_index_t>
+{
+  typedef typename property_map<P, CGAL::vertex_index_t>::type type;
+  typedef typename property_map<P, CGAL::vertex_index_t>::const_type const_type;
+};
+
 } // namespace boost
 
 
@@ -90,6 +97,13 @@ get(boost::vertex_index_t, const Dual<P>& dual)
   return get(CGAL::face_index, dual.primal());
 }
 
+template <typename P>
+typename boost::property_map<P, boost::vertex_index_t>::type
+get(boost::face_index_t, const Dual<P>& dual)
+{
+  return get(CGAL::vertex_index, dual.primal());
+}
+
 
 template <typename P>
 typename boost::graph_traits<CGAL::Dual<P> >::vertices_size_type

Polygon_mesh_processing/examples/Polygon_mesh_processing/connected_component.cpp

@@ -56,6 +56,9 @@ int main(int, char* argv[])
  BOOST_FOREACH(face_descriptor f , faces(sm)){
    std::cout  << f << " in connected component " << fccmap[f] << std::endl;
   }
+ 
+ CGAL::Polygon_mesh_processing::keep_largest_connected_components(sm,2);
 
+ std::cout << "mesh:\n" << sm << std::endl;
   return 0;
 }

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Connected_components.h

@@ -27,7 +27,10 @@
 #include <boost/foreach.hpp>
 #include <boost/graph/filtered_graph.hpp>
 #include <boost/graph/connected_components.hpp>
+#include <boost/property_map/vector_property_map.hpp>
+
 #include <CGAL/boost/graph/iterator.h>
+#include <CGAL/boost/graph/helpers.h>
 
 #include <CGAL/Polyhedron_3.h>
 #include <CGAL/Polyhedron_incremental_builder_3.h>
@@ -337,6 +340,22 @@ mark_connected_components_v2(
 } } // end of namespace internal::corefinement
 
 namespace Polygon_mesh_processing{
+
+  namespace internal {
+    struct LessFirst {
+      typedef std::pair<std::size_t,std::size_t> T;
+      bool operator()(const T& a, const T& b) const {
+        return a.first < b.first;
+      }
+    };
+    struct MoreSecond {
+      typedef std::pair<std::size_t,std::size_t> T;
+      bool operator()(const T& a, const T& b) const {
+        return a.second > b.second;
+      }
+    };
+  } // namespace internal
+
 /*!
  * \ingroup PkgPolygonMeshProcessing
  *  Erases the small connected components and the isolated vertices.
@@ -347,9 +366,126 @@ namespace Polygon_mesh_processing{
 template <class PolygonMesh>
 std::size_t keep_largest_connected_components(PolygonMesh& pmesh, std::size_t nb_components_to_keep)
 {
-  return pmesh.keep_largest_connected_components(nb_components_to_keep);
-}
+  typedef boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
+  typedef boost::graph_traits<PolygonMesh>::face_descriptor face_descriptor;
+  typedef boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
+  typedef boost::graph_traits<PolygonMesh>::edge_descriptor edge_descriptor;
+  boost::vector_property_map<int, boost::property_map<PolygonMesh, boost::face_index_t>::type> face_cc(get(boost::face_index,pmesh));
   
+  int num = connected_components(pmesh,face_cc);
+  if((num == 1)|| (nb_components_to_keep > num) ){
+    return 0;
+  }
+  boost::vector_property_map<bool, boost::property_map<PolygonMesh, boost::vertex_index_t>::type> keep_vertex(get(boost::vertex_index,pmesh)); 
+  BOOST_FOREACH(vertex_descriptor v, vertices(pmesh)){
+    keep_vertex[v] = false;
+  }
+  typedef std::pair<std::size_t, std::size_t> Component;
+  std::vector<Component> component_size(num);
+
+  for(std::size_t i=0; i < num; i++){
+    component_size[i] = std::make_pair(i,0);
+  }
+  BOOST_FOREACH(face_descriptor f, faces(pmesh)){
+    component_size[face_cc[f]].second++;
+  }
+  for(int i=0; i < component_size.size(); ++i){
+    std::cerr << "component " << i << " has " << component_size[i].second << " faces\n";
+  }
+  // we have to sort the range [0, num) by component size
+  internal::MoreSecond ls;
+  std::sort(component_size.begin(), component_size.end(), ls);
+  for(std::size_t i=0; i < num; i++){
+    component_size[i].second = (i < nb_components_to_keep)?1:0;
+  }
+  internal::LessFirst lsinv;
+  std::sort(component_size.begin(), component_size.end(), lsinv);
+  for(std::size_t i=0; i < num; i++){
+    std::cout << i  << " " << component_size[i].first  << " " << component_size[i].second << std::endl;
+  }
+  BOOST_FOREACH(face_descriptor f, faces(pmesh)){
+    face_cc[f] = component_size[face_cc[f]].second;
+    if(face_cc[f] == 1){
+      std::cerr << "keep " << f << std::endl;
+    }
+  }
+  // Now face_cc[f] == 1 means that we want to keep the face
+
+  BOOST_FOREACH(face_descriptor f, faces(pmesh)){
+    if(face_cc[f] == 1){
+      BOOST_FOREACH(halfedge_descriptor h, halfedges_around_face(halfedge(f,pmesh),pmesh)){
+          vertex_descriptor v = target(h,pmesh);
+          keep_vertex[v] = true;
+          std::cout << "keep vertex "<< v << std::endl;
+        }
+    }
+  }
+  BOOST_FOREACH(edge_descriptor e, edges(pmesh)){
+    vertex_descriptor v = source(e,pmesh);
+    vertex_descriptor w = target(e,pmesh);
+    halfedge_descriptor h = halfedge(e,pmesh);
+    halfedge_descriptor oh = opposite(h,pmesh);
+    if(! keep_vertex[v] && ! keep_vertex[w]){
+      // don't care about connectivity
+      // As vertices are not kept the faces and vertices will be removed later
+      remove_edge(e,pmesh);
+    } else if( keep_vertex[v] &&  keep_vertex[w]){
+      face_descriptor fh = face(h,pmesh), ofh = face(oh,pmesh);
+      if(is_border(h,pmesh) && is_border(oh,pmesh)){
+        std::cerr << "null_face on both sides of " << e << " is kept\n";
+      } else if( (face_cc[fh] && is_border(oh,pmesh)) ||
+                 (face_cc[ofh] && is_border(h,pmesh)) ||
+                 (face_cc[fh] && face_cc[ofh]) ){
+        // do nothing
+      } else if(face_cc[fh] && ! face_cc[ofh]){
+        set_face(oh, boost::graph_traits<PolygonMesh>::null_face(), pmesh);
+      } else if(! face_cc[fh] && face_cc[ofh]){
+        set_face(h, boost::graph_traits<PolygonMesh>::null_face(), pmesh);
+      } else {
+        // no face kept
+        assert( ! face_cc[fh] && ! face_cc[ofh]);
+        // vertices pointing to e must change their halfedge
+        if(halfedge(v,pmesh) == oh){
+          set_halfedge(v,prev(h,pmesh),pmesh);
+        }
+        if(halfedge(w,pmesh) == h){
+          set_halfedge(w,prev(oh,pmesh),pmesh);
+        }
+        // shortcut the next pointers as e will be removed
+        set_next(prev(h,pmesh), next(oh,pmesh),pmesh);
+        set_next(prev(oh,pmesh), next(h,pmesh),pmesh);
+        remove_edge(e,pmesh);
+      }
+    } else if( keep_vertex[v] ){
+      if(halfedge(v,pmesh) == oh){
+        set_halfedge(v,prev(h,pmesh),pmesh);
+      }
+      set_next(prev(h,pmesh), next(oh,pmesh),pmesh);
+        remove_edge(e,pmesh);
+    } else {
+      assert (keep_vertex[w]);
+      if(halfedge(w,pmesh) == h){
+        set_halfedge(w,prev(oh,pmesh),pmesh);
+      }
+      set_next(prev(oh,pmesh), next(h,pmesh),pmesh);
+        remove_edge(e,pmesh);
+    }
+  }
+  // We now can remove all vertices and faces not marked as kept
+  BOOST_FOREACH(face_descriptor f, faces(pmesh)){
+    if(face_cc[f] != 1){
+      remove_face(f,pmesh);
+    }
+  }
+  BOOST_FOREACH(vertex_descriptor v, vertices(pmesh)){
+    if(! keep_vertex[v]){
+      remove_vertex(v,pmesh);
+    }
+  }
+
+  return num - nb_components_to_keep;
+}
+
 /*!
  * \ingroup PkgPolygonMeshProcessing
  *  Discovers all the faces in the same connected component as `seed_face` and puts them in `out`.
@@ -404,7 +540,7 @@ struct No_constraint {
 
   template <typename T>
   bool operator[](const T & ) const {
-    return true;
+    return false;
   }
 
   G* g;