Incremental_delaunay.h

// SPDX-License-Identifier: GPL-3.0-or-later
// As an exception, if you have a commercial license for the CGAL package Triangulation,
// you may use this file under the same conditions.
//
// Author(s)    : Michael Kerber, Marc Glisse, (Samuel Hornus for the original CGAL code)
 
#include <CGAL/Delaunay_triangulation.h>
#include <iterator>
#include <vector>
 
namespace Gudhi {
 
template<typename K, typename SimplicialComplex, typename PointRange>
void construct_incremental_delaunay(K const&k, SimplicialComplex& complex, PointRange const& points) {
  using TDS = CGAL::Triangulation_data_structure<
    typename K::Dimension,
    CGAL::Triangulation_vertex<K, typename SimplicialComplex::Vertex_handle>,
    CGAL::Triangulation_ds_full_cell<void, CGAL::TDS_full_cell_mirror_storage_policy> >;
  using Triangulation = CGAL::Delaunay_triangulation<K, TDS>;
  auto pt_it = std::begin(points);
  auto pt_end = std::end(points);
  if(pt_it == pt_end) return;
  int d = k.point_dimension_d_object()(*pt_it);
  Triangulation tri(d);
 
  // local variables pulled out for recycling
  typedef std::vector<typename Triangulation::Full_cell_handle> Full_cell_h_vector;
  Full_cell_h_vector cs;
  std::vector<typename SimplicialComplex::Vertex_handle> splx;
 
  for(int idx = 0; pt_it != pt_end; ++pt_it, ++idx) {
    // Adapted from Delaunay_triangulation::insert(Point)
    auto&& p = *pt_it;
    typename Triangulation::Locate_type lt;
    typename Triangulation::Face f(d);
    typename Triangulation::Facet ft;
    typename Triangulation::Full_cell_handle s = tri.locate(p, lt, f, ft);
    // Adapted from Delaunay_triangulation::insert(Point,Locate_type,Face,Facet,Full_cell_handle)
    if (lt == Triangulation::OUTSIDE_AFFINE_HULL) {
      // nothing disappears
      tri.insert_outside_affine_hull(p)->data() = idx;
    }
    else if (lt == Triangulation::ON_VERTEX) {
      // nothing disappears
    }
    else {
      if(tri.current_dimension() == 1) {
        if(lt == Triangulation::OUTSIDE_CONVEX_HULL) {
          // only an infinite face disappears
        } else {
          complex.insert_simplex_and_subfaces({s->vertex(0)->data(), s->vertex(1)->data(), idx});
        }
        typename Triangulation::Vertex_handle v = tri.tds().insert_in_full_cell(s);
        v->set_point(p);
        v->data() = idx;
      } else { // main case
        // Adapted from Delaunay_triangulation::insert_in_conflicting_cell(Point,Full_cell_handle)
        /* Full_cell_h_vector cs; */ cs.clear();
        std::back_insert_iterator<Full_cell_h_vector> out(cs);
        typename Triangulation::Facet ft = tri.compute_conflict_zone(p, s, out);
        for(auto c : cs) {
          /* std::vector<int> splx; */ splx.clear();
          bool interrupted = false;
          for(int i = 0; i <= tri.current_dimension(); ++i) {
            // We could skip the infinite vertex and output a (new) simplex, but that's not so useful
            if(tri.is_infinite(c->vertex(i))) { interrupted = true; break; }
            splx.push_back(c->vertex(i)->data());
          }
          if (interrupted) continue;
          // We could put idx at the beginning of the vector once instead of reinserting it each time
          splx.push_back(idx);
          complex.insert_simplex_and_subfaces(splx);
        }
        tri.insert_in_hole(p, cs.begin(), cs.end(), ft)->data() = idx;
      }
    }
  }
  // We have added all the conflict cells, but we are missing some current cells.
  // There will be a lot of redundancy though.
  for(auto it = tri.finite_full_cells_begin(); it != tri.finite_full_cells_end(); ++it) {
    /* std::vector<int> splx; */ splx.clear();
    for(int i = 0; i <= tri.current_dimension(); i++) {
      splx.push_back(it->vertex(i)->data());
    }
    complex.insert_simplex_and_subfaces(splx);
    // Or we could call insert_simplex_and_subfaces on a transform_iterator.
  }
}
 
/* Strategies
 1. Current
 - every time a finite simplex disappears, output the conflict simplex
 - output the final triangulation at the end (a bit costly)
 2. Alternative from Kerber's function_delaunay
 - every time a full-dimensional simplex disappears, output the finite faces of the conflict simplex
 - this only works with some assumptions, for instance that the first d+1 points are affinely independent and
   that the simplex they define is not in the final Delaunay triangulation.
 3. Mix?
 - we could use 1. until the triangulation reaches the ambient dimension then switch to 2.
 - we could pretend the current dimension is going to be the final one, use 2, but fix things up every time
   we insert a point outside the affine hull. Probably not very interesting, since triangulations of non-full
   dimension should be rare.
*/
 
/* Possibilities:
 - output conflict cells, new boundary cells, and final cells separately in case a user wants to handle them differently?
 - we could remove duplicates in `points` so the numbering becomes contiguous, but then we have to output the new `points`.
 - use less undocumented functions. Calling `insert` shouldn't cost that much more.
*/
} // namespace Gudhi