Lloyd_optimize_mesh, constrained delaunay triangulation and assertion violation

[yg42]

Issue Details

Dear all, I submit my problem to your wisdom...

When running the next code, I get an assertion violation:

terminate called after throwing an instance of 'CGAL::Assertion_exception' what(): CGAL ERROR: assertion violation! Expr: v == V[2] File: /home/yann/Documents/Compilation/cgal/TDS_2/include/CGAL/Triangulation_ds_face_base_2.h Line: 167

Is there a bug or is it a misunderstanding or misuse on my part? in this case, I am really sorry to bother. Thank you very much for your time.

Source Code

Here is the source code I tried to minimize.

//#define CGAL_MESH_2_OPTIMIZER_VERBOSE
//#define CGAL_MESH_2_OPTIMIZERS_DEBUG
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/Delaunay_mesher_2.h>
#include <CGAL/Delaunay_mesh_face_base_2.h>
#include <CGAL/Delaunay_mesh_vertex_base_2.h>

#include <CGAL/lloyd_optimize_mesh_2.h>
#include <CGAL/point_generators_2.h>
#include <CGAL/random_selection.h>

#include <limits.h>

#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/Delaunay_triangulation_adaptation_traits_2.h>
#include <CGAL/Delaunay_triangulation_adaptation_policies_2.h>

#include <iostream>

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Delaunay_mesh_vertex_base_2<K>                Vb;
typedef CGAL::Delaunay_mesh_face_base_2<K>                  Fb;
typedef CGAL::Triangulation_data_structure_2<Vb, Fb>        Tds;
typedef CGAL::Constrained_Delaunay_triangulation_2<K, Tds>  CDT;
typedef CGAL::Delaunay_triangulation_2<K>                   DT;

typedef CDT::Point Point;
typedef CGAL::Creator_uniform_2<double,Point>            Creator;

typedef CDT::Vertex_handle Vertex_handle;

#define RADIUS 5.5
#define RADIUS_BIRTH 5

void cell_death(CDT &cdt)
{
  CDT::Finite_vertices_iterator it;
  // CDT::Vertex_handle v; // equivalent à l'iterateur

  std::vector<CDT::Vertex_handle> vector_vertices;
  std::vector<CDT::Vertex_handle> dead_cells;
  std::vector<CDT::Vertex_handle>::iterator it_cells;
  int count=0;
  for (it = cdt.finite_vertices_begin();
       it!= cdt.finite_vertices_end();
       ++it)
    {

      // on ne garde que les points à l'intérieur des contraintes
      // delete only points without constraints
      if (!cdt.are_there_incident_constraints(it))
        {
          vector_vertices.push_back(it);
          count++;
        }

    }

  int n =10;

  random_selection(vector_vertices.begin(), vector_vertices.end(), n,
           std::back_inserter(dead_cells));

  for (it_cells = dead_cells.begin(); it_cells!= dead_cells.end(); ++it_cells)
    cdt.remove(*it_cells);

}

void cell_birth(CDT &cdt)
{
  // insert random points
  CGAL::Random_points_in_disc_2<Point,Creator> g(RADIUS_BIRTH);
  int new_cells = 10;
  for (int i=0; i<new_cells; ++i)
    {
      cdt.insert(*g++);
    }

}

int startSimu(int nb_points, int nb_iter)
{

  CDT cdt;
  // insert random points at first
  CGAL::Random_points_in_disc_2<Point,Creator> g(RADIUS);
  for (int i=0; i<nb_points; ++i)
    {
      cdt.insert(*g++);
    }

  // loop on birth death and lloyd optimization
    for (int n=0; n<100; ++n)
    {

      cell_birth(cdt);
      cell_death(cdt);
      lloyd_optimize_mesh_2(cdt,
                    CGAL::parameters::max_iteration_number = nb_iter,
                    CGAL::parameters::mark = false);
    }

  return 0;
}

/// main function
int main(int argc, char**argv)
{
  // arguments
  int nb_points=3000; // nb of simulated points
  int nb_iter=5; // nb of iterations

  return startSimu(nb_points, nb_iter);
}

Environment

• Operating system (Windows/Mac/Linux, 32/64 bits): Linux ubuntu 20.10, 64 bits
• Compiler: gcc (Ubuntu 10.2.0-13ubuntu1) 10.2.0
• Release or debug mode: both
• Specific flags used (if any):
• CGAL version: from git 5.2.x branch
• Boost version: 1.71.0-6ubuntu9
• Other libraries versions if used (Eigen, TBB, etc.):

cmake_minimum_required(VERSION 3.1...3.14)
project( life )
set(CMAKE_CXX_STANDARD 14)
find_package(CGAL REQUIRED QUIET OPTIONAL_COMPONENTS Core Qt5)

set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall" )

if(CGAL_FOUND AND CGAL_Qt5_FOUND)
  #required to use basic_viewer
  add_definitions(-DCGAL_USE_BASIC_VIEWER -DQT_NO_KEYWORDS)

  create_single_source_cgal_program( "mwe.cpp" )

else()
  message("error with Qt5")
endif()

To reproduce:

mkdir debug
cd debug
cmake -g .. -DCMAKE_BUILD_TYPE=Debug -DCGAL_DIR=~/Documents/Compilation/cgal
./mwe

[janetournois]

Hello @yg42 ,

thank you for your report, and sorry for the late reply.

This issue comes from the fact that random_selection() can collect twice (or more) the same vertex in dead_cells. So the line cdt.remove(*it_cells); is likely to fail at removing a vertex that has already been removed earlier in the same loop. You can use a std::set to store the "dead vertices" for example to make sure that your vertex set is unique.

[yg42]

Thank you very much for your time. I can thus propose to add in Generator/include/CGAL/random_selection.h the following code. I am really sorry I do not have time for proposing a pull request.

template <class RandomAccessIterator, class Size, class OutputIterator,
          class Random>
OutputIterator unique_random_selection( RandomAccessIterator first,
                    RandomAccessIterator last,
                    Size n,
                    OutputIterator result,
                    Random& rnd)
   // choose n random (unique) values from the range [`first',`last')
   // and write it to `result'. Each item has an equal probability to
   // be chosen. In case n is greater than range, this is silently ignored
{
  std::ptrdiff_t m = last - first;
  std::vector<int> numbers;
  for (int i=0; i<m; i++) {
    numbers.push_back(i);
  }
  std::iota(numbers.begin(), numbers.end(), 0);

  // ignore in case to many items are requested
  if (n > m)
    n = m;

  // random shuffle, then get only n first items.
  CGAL::cpp98::random_shuffle(numbers.begin(), numbers.end(), rnd);
  for ( Size i = 0; i < n; i++) {
    *result++ = first[ numbers[i]];
  }
  return result;
}

template <class RandomAccessIterator, class Size, class OutputIterator>
OutputIterator unique_random_selection( RandomAccessIterator first,
                    RandomAccessIterator last,
                    Size n,
                    OutputIterator result)
{
    return unique_random_selection( first, last, n, result, get_default_random());
}