stitch_cycles.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>

#include <CGAL/Polygon_mesh_processing/stitch_borders.h>
#include <CGAL/boost/graph/copy_face_graph.h>

#include <iostream>
#include <fstream>

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Surface_mesh<K::Point_3> Mesh;

template <class PolygonMesh, class OutputIterator>
OutputIterator
extract_boundary_cycles(const PolygonMesh& pm, OutputIterator out)
{
  typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
  std::vector<halfedge_descriptor> border_hedges;

  BOOST_FOREACH(halfedge_descriptor h, halfedges(pm))
  {
    if (is_border(h, pm))
      border_hedges.push_back(h);
  }

  // use a dynamic property but mesh is const. Request index map?
  boost::unordered_set<halfedge_descriptor> handled;
  BOOST_FOREACH(halfedge_descriptor h, border_hedges)
  {
    if (!handled.insert(h).second) continue;
    *out++=h;
    halfedge_descriptor hn=next(h, pm);
    while(hn!=h)
    {
      handled.insert(hn);
      hn=next(hn, pm);
    }
  }
  return out;
}


int main(int, char** argv)
{
  Mesh sm;
  std::ifstream(argv[1]) >> sm;

  std::vector<Mesh::Halfedge_index> boundary_cycles;
  extract_boundary_cycles(sm, std::back_inserter(boundary_cycles));

  std::cout << "Found " << boundary_cycles.size() << " boundary cycles\n";
  
  BOOST_FOREACH(Mesh::Halfedge_index h, boundary_cycles)
  {
    Mesh::Halfedge_index hn=next(h, sm);
    bool do_stitching=false;
    while(hn!=h)
    {
      if ( sm.point( source(hn, sm) ) == sm.point( target(next(hn, sm), sm) ) )
      {
        do_stitching=true;
        break;
      }
      hn=next(hn, sm);
    }

    if (!do_stitching) continue;

    std::vector< std::pair<Mesh::Halfedge_index, Mesh::Halfedge_index> > hedges_to_stitch;
    
    Mesh::Halfedge_index hnn=next(hn, sm);
    do{
      hedges_to_stitch.push_back( std::make_pair(hn, hnn) );

      if (next(hnn, sm)==hn) break;

      hn=prev(hn, sm);
      hnn = next(hnn, sm);
      
      if ( sm.point( source(hn, sm) ) != sm.point( target(hnn, sm) ) )
        do_stitching = false;
      
    }while(do_stitching);

    if (!do_stitching) continue;

    std::cout << "stitching a cycle\n";
    
    CGAL::Polygon_mesh_processing::stitch_borders(sm, hedges_to_stitch);

    std::ofstream("out.off") << std::setprecision(17) << sm;
  }


  return 0;
}