Secure usage of Cell_handle in a CGAL Delaunay triangulation after point insertion

Question

I'm planning to write an algorithm that will use CGAL Delaunay triangulation data structure. Basically I need to insert some point into the triangulation, save reference to some cells, and then make some other insertion.

I'm wondering how can I store reference to cells that are not invalidated after insertion of new points in triangulation?

It's seems to me that Cell_handle is just a pointer to an internal structure, so it's dangerous to store it due to reallocation of internal container. In the other hand I can see no way in Triangulation_3 interface to store an index from a Cell_handle.

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Triangulation_vertex_base_3<K>             Vb;
typedef CGAL::Triangulation_hierarchy_vertex_base_3<Vb>  Vbh;
typedef CGAL::Triangulation_data_structure_3<Vbh>        Tds;
typedef CGAL::Delaunay_triangulation_3<K,Tds>            Dt;
typedef CGAL::Triangulation_hierarchy_3<Dt>              Dh;
typedef Dh::Vertex_iterator Vertex_iterator;
typedef Dh::Vertex_handle   Vertex_handle;
typedef Dh::Point           Point;

int main(){

   Dh T;
   for(int i = 0; i < 100; ++i)
      T.insert(Point(rand()%1000,rand()%1000,rand()%1000));

   assert( T.is_valid() );
   assert( T.number_of_vertices() == 100 );
   assert( T.dimension() == 3 );

   typedef Dh::Cell_iterator CellIterator;

   std::vector<Dh::Cell_handle> hnd;
   CellIterator itEnd = T.finite_cells_end();
   for(CellIterator it = T.finite_cells_begin(); it!=itEnd; ++it){
      const int dist = std::distance(T.cells_begin(),it);
      hnd.push_back(it);
   }

   const int newP(1000);
   for(int i = 0; i < newP; ++i)
      T.insert(Point(rand()%1000,rand()%1000,rand()%1000));

   int finiteC(0),infiniteC(0);
   for(int i = 0; i < hnd.size(); ++i){
      const int dist = std::distance(T.cells_begin(),hnd[i]);
      if(T.is_infinite(hnd[i]))
     ++infiniteC;
  else
     ++finiteC;
   }
   assert( T.is_valid() );
   return 0;
}

This code systematically crash but, and this is really strange to me, if I change newP to 10000, this code magically works.

Can someone explain me how to handle this problem?

Answer 1

Since a cell can disappear during insertion of a new point, the handle you have saved are not guarantee to point on what you expect.

You have a crash because you use the triangulation hierarchy that internally creates and remove cells in the internal container. If you use CGAL::Delaunay_triangulation_3, you will not have the crash.

For your problem, you should store a quadruplet of Vertex_handleS and use the is_cell function (documented here).