Retrieve vertices from CGAL's Delaunay Constrained Triangulation

Question

I have a list of 2D points, which describe a simple polygon. I compute the Delaunay constrained thanks to CGAL.

    Polygon_2 polygon;
    //Create a polygon from a vector
    for(int j=0; j < (*itVectorPoint2D).size(); j++) {
        polygon.push_back(Point((*itVectorPoint2D)[j].x,(*itVectorPoint2D)[j].y));
    }

    //Insert the polygon into a constrained triangulation
    CDT cdt;
    insert_polygon(cdt,polygon);
    //Mark facets that are inside the domain bounded by the polygon
    mark_domains(cdt);

    for (CDT::Finite_faces_iterator fit=cdt.finite_faces_begin(); fit!=cdt.finite_faces_end();++fit)
    {
        if (fit->info().in_domain()) {
            std::cout << (*(fit->vertex(0))) << ";" << (*(fit->vertex(1))) << ";" << (*(fit->vertex(2))) << std::endl;
        }

    }

My problem is that I want to retrieve every triangle with the triplet of indices from the initial vector.For example, for a square define by 0:(0,0) 1:(1,0) 2:(1,1) and 3:(0,1), I want to retrieve (0, 1, 2) and (0, 2, 3). I could retrieve the coordinates of the current vertex and compare them to every point in the vector, but it's ugly I think (and inefficient, possible source of problem ?).

Is there a way to do so ?

Answer 1

The following works fine with CGAL 4.2. Note that with CGAL 4.3 there will be a dedicated function doing the job in a more integrated way like it is done here for the Delaunay_triangulation_2 class.

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/Triangulation_vertex_base_with_info_2.h>
#include <CGAL/Spatial_sort_traits_adapter_2.h>
#include <vector>

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;

typedef CGAL::Triangulation_vertex_base_with_info_2<unsigned, K> Vb;
typedef CGAL::Constrained_triangulation_face_base_2<K>           Fb;
typedef CGAL::Triangulation_data_structure_2<Vb,Fb>              TDS;
typedef CGAL::Exact_predicates_tag                               Itag;
typedef CGAL::Constrained_Delaunay_triangulation_2<K, TDS, Itag> CDT;
typedef CDT::Point          Point;
typedef CGAL::Spatial_sort_traits_adapter_2<K,Point*> Search_traits;

template <class InputIterator>
void insert_with_info(CDT& cdt, InputIterator first,InputIterator last)
{
  std::vector<std::ptrdiff_t> indices;
  std::vector<Point> points;
  std::ptrdiff_t index=0;

  for (InputIterator it=first;it!=last;++it){
    points.push_back( *it);
    indices.push_back(index++);
  }

  CGAL::spatial_sort(indices.begin(),indices.end(),Search_traits(&(points[0]),cdt.geom_traits()));

  CDT::Vertex_handle v_hint;
  CDT::Face_handle hint;
  for (typename std::vector<std::ptrdiff_t>::const_iterator
    it = indices.begin(), end = indices.end();
    it != end; ++it){
    v_hint = cdt.insert(points[*it], hint);
    if (v_hint!=CDT::Vertex_handle()){
      v_hint->info()=*it;
      hint=v_hint->face();
    }
  }
}

int main()
{
  std::vector< Point > points;
  points.push_back( Point(0,0) );
  points.push_back( Point(1,0) );
  points.push_back( Point(0,1) );
  points.push_back( Point(14,4) );
  points.push_back( Point(2,2) );
  points.push_back( Point(-4,0) );


  CDT cdt;
  insert_with_info(cdt, points.begin(), points.end());

  CGAL_assertion( cdt.number_of_vertices() == 6 );

  // check that the info was correctly set.
  CDT::Finite_vertices_iterator vit;
  for (vit = cdt.finite_vertices_begin(); vit != cdt.finite_vertices_end(); ++vit)
    if( points[ vit->info() ] != vit->point() ){
      std::cerr << "Error different info" << std::endl;
      exit(EXIT_FAILURE);
    }
  std::cout << "OK" << std::endl;
}