terrain2sm.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/Projection_traits_xy_3.h>
#include <CGAL/Triangulation_face_base_with_info_2.h>
#include <CGAL/Polygon_2.h>
#include <CGAL/Surface_mesh.h>

// needed for copy_face_graph
#include <CGAL/boost/graph/Face_filtered_graph.h>
#include <CGAL/boost/graph/graph_traits_Constrained_delaunay_triangulation_2.h>
#include <boost/property_map/function_property_map.hpp>

// needed for mesh simplification
#include <CGAL/Surface_mesh_simplification/edge_collapse.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Count_ratio_stop_predicate.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Bounded_normal_change_placement.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_policies.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Constrained_placement.h>

#include <iostream>
#include <fstream>


struct FaceInfo2
{
  FaceInfo2(){}
  int nesting_level;

  bool in_domain(){
    return nesting_level%2 == 1;
  }
};


typedef CGAL::Exact_predicates_inexact_constructions_kernel       K;
typedef CGAL::Projection_traits_xy_3<K>  Gt;
typedef CGAL::Triangulation_vertex_base_2<Gt>                     Vb;
typedef CGAL::Triangulation_face_base_with_info_2<FaceInfo2,Gt>   Fbb;
typedef CGAL::Constrained_triangulation_face_base_2<Gt,Fbb>       Fb;
typedef CGAL::Triangulation_data_structure_2<Vb,Fb>               TDS;
typedef CGAL::Exact_predicates_tag                                Itag;
typedef CGAL::Constrained_Delaunay_triangulation_2<Gt, TDS, Itag> CDT;
typedef CDT::Point                                                Point;
typedef CGAL::Polygon_2<Gt>                                       Polygon_2;
typedef CDT::Face_handle                                          Face_handle;

typedef CGAL::Face_filtered_graph<CDT>                            Domain_face_graph;
typedef CGAL::Surface_mesh<Point>                                 Surface_mesh;
typedef boost::graph_traits<Surface_mesh>::halfedge_descriptor    halfedge_descriptor;
typedef boost::graph_traits<Surface_mesh>::edge_descriptor        edge_descriptor;


namespace SMS = CGAL::Surface_mesh_simplification;

// Garland&Heckbert simplification policies
typedef typename SMS::GarlandHeckbert_policies<Surface_mesh, K>   GH_policies;
typedef typename GH_policies::Get_cost                            GH_cost;
typedef typename GH_policies::Get_placement                       GH_placement;
typedef SMS::Bounded_normal_change_placement<GH_placement>        Bounded_GH_placement;

struct Border_is_constrained_edge_map
{
  const Surface_mesh* sm_ptr;
  typedef edge_descriptor                                       key_type;
  typedef bool                                                  value_type;
  typedef value_type                                            reference;
  typedef boost::readable_property_map_tag                      category;

  Border_is_constrained_edge_map(const Surface_mesh& sm) : sm_ptr(&sm) {}

  friend bool get(Border_is_constrained_edge_map m, const key_type& edge) {
    return CGAL::is_border(edge, *m.sm_ptr);
  }
};

typedef SMS::Constrained_placement<GH_placement,
                                   Border_is_constrained_edge_map > Placement;


void
mark_domains(CDT& ct,
             Face_handle start,
             int index,
             std::list<CDT::Edge>& border )
{
  if(start->info().nesting_level != -1){
    return;
  }
  std::list<Face_handle> queue;
  queue.push_back(start);

  while(! queue.empty()){
    Face_handle fh = queue.front();
    queue.pop_front();
    if(fh->info().nesting_level == -1){
      fh->info().nesting_level = index;
      for(int i = 0; i < 3; i++){
        CDT::Edge e(fh,i);
        Face_handle n = fh->neighbor(i);
        if(n->info().nesting_level == -1){
          if(ct.is_constrained(e)) border.push_back(e);
          else queue.push_back(n);
        }
      }
    }
  }
}

//explore set of facets connected with non constrained edges,
//and attribute to each such set a nesting level.
//We start from facets incident to the infinite vertex, with a nesting
//level of 0. Then we recursively consider the non-explored facets incident
//to constrained edges bounding the former set and increase the nesting level by 1.
//Facets in the domain are those with an odd nesting level.
void
mark_domains(CDT& cdt)
{
  for(CDT::Face_handle f : cdt.all_face_handles()){
    f->info().nesting_level = -1;
  }

  std::list<CDT::Edge> border;
  mark_domains(cdt, cdt.infinite_face(), 0, border);
  while(! border.empty()){
    CDT::Edge e = border.front();
    border.pop_front();
    Face_handle n = e.first->neighbor(e.second);
    if(n->info().nesting_level == -1){
      mark_domains(cdt, n, e.first->info().nesting_level+1, border);
    }
  }
}


int main( )
{
  //construct two non-intersecting nested polygons
  Polygon_2 polygon1;
  polygon1.push_back(Point(0,0,0));
  polygon1.push_back(Point(2,0,0));
  polygon1.push_back(Point(2,2,0));
  polygon1.push_back(Point(0,2,0));
  Polygon_2 polygon2;
  polygon2.push_back(Point(0.5,0.5,0));
  polygon2.push_back(Point(1.5,0.5,0));
  polygon2.push_back(Point(1.5,1.5,0));
  polygon2.push_back(Point(0.5,1.5,0));

  Polygon_2 polygon3;
  polygon3.push_back(Point(10,0,0));
  polygon3.push_back(Point(12,0,0));
  polygon3.push_back(Point(12,2,0));
  polygon3.push_back(Point(10,2,0));

  //Insert the polygons into a constrained triangulation
  CDT cdt;
  cdt.insert_constraint(polygon1.vertices_begin(), polygon1.vertices_end(), true);
  cdt.insert_constraint(polygon2.vertices_begin(), polygon2.vertices_end(), true);
  cdt.insert_constraint(polygon3.vertices_begin(), polygon3.vertices_end(), true);

  cdt.insert(Point(-1, -1, 0)); // outside the domain
  cdt.insert(Point(1, 1, 0)); // outside the domain
  cdt.insert(Point(0.25, 0.25, 0.2));
  cdt.insert(Point(11, 1, 0.3));

  //Mark facets that are inside the domain bounded by the polygon
  mark_domains(cdt);

  Surface_mesh sm;

  Domain_face_graph dfg(cdt, true, boost::make_function_property_map<Face_handle, bool>([](Face_handle f ){ return f->info().in_domain(); }));
  CGAL::copy_face_graph(dfg, sm);

  SMS::Count_stop_predicate<Surface_mesh> stop(0); // Contract the surface mesh as much as possible
  Border_is_constrained_edge_map bem(sm);

  GH_policies gh_policies(sm);
  const GH_cost& gh_cost = gh_policies.get_cost();
  const GH_placement& gh_placement = gh_policies.get_placement();
  Bounded_GH_placement bgh_placement(gh_placement);


  SMS::edge_collapse(sm, stop,
                     CGAL::parameters::edge_is_constrained_map(bem)
                     .get_placement(Placement(bem, bgh_placement)));

  std::ofstream out("out.off");
  out << sm << std::endl;

  return 0;
}