How do I iterate over faces in CGAL

Question

I am trying to use CGAL to do some Delaunay triangulation. I used one of the CGAL samples to compute a triangulation which includes a height field attribute.

The problem I have having is that I have no idea how to get the resulting triangulation. I figured out how to get the face_iterator, but I don't know what to do from there. What I'm hoping to get is an index into the point array for each of the 3 points on each triangle.

I'm having trouble wading through all of the nested templates:

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Triangulation_euclidean_traits_xy_3.h>
#include <CGAL/Delaunay_triangulation_2.h>

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Triangulation_euclidean_traits_xy_3<K> Gt;
typedef CGAL::Delaunay_triangulation_2<Gt> Delaunay;
typedef K::Point_3 Point;

int main()
{
    //initialize the points with some trivial data
    std::vector<Point> pts;
    pts.push_back(Point(1., 2., 3.));
    pts.push_back(Point(2., 2., 3.));
    pts.push_back(Point(1., 3., 3.));
    pts.push_back(Point(4., 2., 3.));    

    //create a delaunay triangulation
    Delaunay dt;
    dt.insert(pts.begin(), pts.end());

    //iterate through the faces
    Delaunay::Finite_faces_iterator it;
    for (it = dt.finite_faces_begin(); it != dt.finite_faces_end(); it++)
    {
        //What do I do here??
    }

    return 0;
}

Heuser

it.something ... does code completion help at all? – Thorax

Visual Studio 2008 is unable to perform code complete on this code :( – Heuser

Highlight Delaunay::Finite_faces_iterator and go to definition. Trace the .h file and see what functions are available. It will be it.<something> – Thorax

a first problem with your code is that you use Delaunay_triangulation_2<Gt> instead of Delaunay_triangulation_3<Gt> – Euonymus

...or use 2D points if you use Delaunay_triangulation_2<Gt>... – Loriannlorianna

CGAL has confusing API and samples. If not for this question, I would've been stuck on this for hours. – Berkey

Answer 1

You can use Delaunay::triangle to convert from a face (iterator) to the corresponding triangle. This is tested under CGAL 3.8:

// points.cin contains point pairs, e.g.,
// 3 5 
// 0 0
// 1 9
// ...
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <fstream>

typedef CGAL::Exact_predicates_exact_constructions_kernel K;
typedef CGAL::Delaunay_triangulation_2<K> Delaunay;
typedef K::Point_2   Point;

int main()
{
  std::ifstream in("points.cin");
  std::istream_iterator<Point> begin(in);
  std::istream_iterator<Point> end;

  Delaunay dt;
  dt.insert(begin, end);

  Delaunay::Finite_faces_iterator it;
  for (it = dt.finite_faces_begin(); it != dt.finite_faces_end(); it++)
  {
    std::cout << dt.triangle(it) << std::endl;
  }

  return 0;
}

Complement

Answer 2

A vertex of a triangle can be accessed using dt.triangle(it)[idx], where it is a faces iterator and idx is vertex number (0,1 or 2). In the example below, a vertex is a Point_2 object, its cartesian coordinates can be accessed using x() and y() methods.

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Triangulation_euclidean_traits_2.h>
#include <CGAL/Delaunay_triangulation_2.h>

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Triangulation_euclidean_traits_2<K> Gt;
typedef CGAL::Delaunay_triangulation_2<Gt> Delaunay;

typedef K::Point_2 Point_2;
typedef std::vector<Point_2> Points;

int main()
{
    Points points;
    points.push_back(Point_2(0,0));
    points.push_back(Point_2(0,7));
    points.push_back(Point_2(7,0));
    points.push_back(Point_2(7,7));

    Delaunay dt(points.begin(), points.end());

    // Print Cartesian coordinates of vertices of triangles in 2D Delaunay triangulation
    for (Delaunay::Finite_faces_iterator it = dt.finite_faces_begin(); it != dt.finite_faces_end(); it++)
    {
        std::cout << " " << dt.triangle(it)[0].x() << " " << dt.triangle(it)[0].y() << " ";
        std::cout << " " << dt.triangle(it)[1].x() << " " << dt.triangle(it)[1].y() << " ";
        std::cout << " " << dt.triangle(it)[2].x() << " " << dt.triangle(it)[2].y() << " ";
        std::cout << std::endl << "-------------------" << std::endl;
    }
    return 0;
}

Faithless

Answer 3

Here is an example from Google. Finite_faces_iterator was typedefed.

  Interval_skip_list isl;
  for(Finite_faces_iterator fh = dt.finite_faces_begin();
      fh != dt.finite_faces_end();
      ++fh){
    isl.insert(Interval(fh));
  }
  std::list<Interval> level;
  isl.find_intervals(50, std::back_inserter(level));
  for(std::list<Interval>::iterator it = level.begin();
      it != level.end();
      ++it){
    std::cout << dt.triangle(it->face_handle()) << std::endl;
  }

This does not do what you want, but gives you an example of what can be done with an iterator.

Thorax

The iterators themselves are stl iterators. The issue is that I can't figure out what the iterator is pointing to (presumably a face or facet or something) – Heuser

Answer 4

If you want a really extended example of how to do exactly what you want, take a look at the source to the X-Plane scenery tools from here: http://scenery.x-plane.com/code.php

By extended example, I mean a couple of hundred thousand lines, but there's uses of nearly everything CGAL can do with Delaunay triangulations and extended attributes in there.

Snooze

Answer 5

Well I just had a similar problem where I did a lot of research (mainly because I did not have any knowledge of C++). I wanted to be able to print triangles by it's vertex integer representation. Here it is how it looks:

#include <CGAL/Surface_mesh_default_triangulation_3.h>
#include <CGAL/Complex_2_in_triangulation_3.h>
#include <CGAL/make_surface_mesh.h>
#include <CGAL/Implicit_surface_3.h>

// This is the file where you can look for an example of iterating, geting basic vertex positions, outputing triangles
// #include <CGAL/IO/Complex_2_in_triangulation_3_file_writer.h>

// default triangulation for Surface_mesher
typedef CGAL::Surface_mesh_default_triangulation_3 Tr;

// c2t3
typedef CGAL::Complex_2_in_triangulation_3<Tr> C2t3;

typedef Tr::Geom_traits GT;
typedef GT::Sphere_3 Sphere_3;
typedef GT::Point_3 Point_3;
typedef GT::FT FT;

typedef FT (*Function)(Point_3);

typedef CGAL::Implicit_surface_3<GT, Function> Surface_3;

// This already have been defined
//typedef typename C2t3::Triangulation Tr;
typedef typename Tr::Vertex_handle Vertex_handle;
typedef typename Tr::Finite_vertices_iterator Finite_vertices_iterator;
typedef typename Tr::Finite_facets_iterator Finite_facets_iterator;

typedef typename Tr::Point Point;


FT sphere_function (Point_3 p) {
  const FT x = p.x();
  const FT y = p.y();
  const FT z = p.z();

  //const FT x2=p.x()*p.x(), y2=p.y()*p.y(), z2=p.z()*p.z();
  const FT a = 2;
  const FT b = 1;
  const FT c = 1.5;
  return x*x/a/a + y*y/b/b + z*z/c/c -1;
}

int main() {
  Tr tr;            // 3D-Delaunay triangulation
  C2t3 c2t3 (tr);   // 2D-complex in 3D-Delaunay triangulation

  // defining the surface
  Surface_3 surface(sphere_function,             // pointer to function
                    Sphere_3(CGAL::ORIGIN, 2.)); // bounding sphere
  // Note that "2." above is the *squared* radius of the bounding sphere!

  // defining meshing criteria
  CGAL::Surface_mesh_default_criteria_3<Tr> criteria(30.,  // angular bound
                                                     0.1,  // radius bound
                                                     0.1); // distance bound
  // meshing surface
  CGAL::make_surface_mesh(c2t3, surface, criteria, CGAL::Non_manifold_tag());

  std::cout << "Final number of points: " << tr.number_of_vertices() << "\n";

  // Here should be the main code

  Tr& tr2 = c2t3.triangulation();

  std::map<Vertex_handle, int> V;
  int inum = 0;
  Finite_vertices_iterator vit = tr2.finite_vertices_begin();
  while(vit != tr2.finite_vertices_end()) {

    // making an integer representation of vertex pointers
    V[vit] = inum++;

    // obtaining vertex positions from vertex pointer vit
    Point p = static_cast<Point>(vit->point());
    std::cout << p.x() << " " << p.y() << " " << p.z() << std::endl;

    ++vit;
  }

  Finite_facets_iterator fit = tr2.finite_facets_begin();

  while (fit != tr2.finite_facets_end()) {

    typename Tr::Cell_handle cell = fit->first;
    const int& index = fit->second;

    int index1 = V[cell->vertex(tr.vertex_triple_index(index, 0))];
    int index2 = V[cell->vertex(tr.vertex_triple_index(index, 1))];
    int index3 = V[cell->vertex(tr.vertex_triple_index(index, 2))];

    std::cout << index1 << " " << index2 << " " << index3 << std::endl;
    ++fit;
  }

}

compile it with (if mesh_implicit_function is source, object file, and executable):

c++   -DCGAL_USE_GMP -DCGAL_USE_MPFR -DCGAL_USE_ZLIB -frounding-math -o mesh_an_implicit_function.cpp.o -c mesh_an_implicit_function.cpp
c++ mesh_an_implicit_function.cpp.o  -o mesh_an_implicit_function -lmpfr -lgmp -lCGAL -lboost_thread

Residuary