How can the output of iterating through a CGAL Constrained_triangulation_plus_2 be made repeatable on different runs?

Question

I’m currently working with CGAL’s 2D constrained triangulation. After inserting polylines into a constrained triangulation, I noticed that, when iterating through them using a Constraint_iterator, the polylines are not always traversed in the same order, even when I run the same code twice on the same machine without recompiling between runs. The standard output of the code below shows the same polylines, with the exact same points, but not always printed in the same order between different runs. Is there a way to have these polylines always iterated through in the same order? Repeatable output, if possible, would help as I am currently trying to debug a larger CGAL project.

I also noticed that this inconsistent behavior happens rarely; sometimes, I have to run the demo about a dozen times before the output shows any inconsistency or run the demo with different input polygons before trying again. Simpler polygons seem to rarely yield this non-deterministic iteration (hence, the complicated polygon example in the following code).

The example I used has polygons with shared boundaries and holes. Here is the image of the example polygons used. I need this polyline traversal to debug code that also calls CGAL’s simplify() function. Therefore, my demo code is modified from the documentation of CGAL’s polyline simplification at https://doc.cgal.org/latest/Polyline_simplification_2/index.html.

#include <iostream>
#include <fstream>
#include <sstream>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polygon_2.h>
#include <CGAL/Polygon_with_holes_2.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/Constrained_triangulation_plus_2.h>
#include <CGAL/IO/WKT.h>

typedef CGAL::Simple_cartesian<double> K;
typedef CGAL::Polygon_2<K> Polygon_2;
typedef CGAL::Polygon_with_holes_2<K> Polygon_with_holes_2;

typedef CGAL::Constrained_Delaunay_triangulation_2<K, CGAL::Default, CGAL::Exact_predicates_tag> CDT;
typedef CGAL::Constrained_triangulation_plus_2<CDT>     CT;
typedef CT::Constraint_iterator             Constraint_iterator;
typedef CT::Points_in_constraint_iterator   Points_in_constraint_iterator;

int main(int argc, char* argv[])
{
  // Uncomment if using an external WKT file
  // std::ifstream ifs( (argc==1)?"data/polygon.wkt":argv[1]);

  // Comment out if using an external WKT file
  std::stringstream ifs("POLYGON ((4.24383 50.8196, 4.38314 50.7637, 4.48228 50.793, 4.43286 50.8947, 4.29449 50.8892, 4.24383 50.8196))\nPOLYGON ((5.68787 50.8119, 5.64166 50.8642, 5.75882 50.9517, 5.75839 51.035, 5.85578 51.1446, 5.24219 51.3052, 5.07117 51.3935, 4.66969 51.4264, 4.53818 51.4824, 4.38484 51.4493, 4.43169 51.375, 4.24205 51.354, 4.07063 51.2507, 3.88635 51.2002, 3.59082 51.3044, 3.44824 51.2415, 3.36613 51.3693, 3.07833 51.3027, 2.5455 51.089, 2.71077 50.8136, 2.86315 50.7083, 3.019 50.7736, 3.17839 50.7561, 3.56821 50.7281, 3.68874 50.7747, 3.8798 50.7515, 3.93145 50.6895, 4.30794 50.7004, 4.46355 50.7551, 4.64354 50.746, 4.74861 50.8074, 4.98168 50.7698, 5.09584 50.7034, 5.38958 50.748, 5.47863 50.7235, 5.68787 50.8119), (4.24383 50.8196, 4.29449 50.8892, 4.43286 50.8947, 4.48228 50.793, 4.38314 50.7637, 4.24383 50.8196))\nPOLYGON ((5.68206 50.7575, 5.68787 50.8119, 5.47863 50.7235, 5.38958 50.748, 5.09584 50.7034, 4.98168 50.7698, 4.74861 50.8074, 4.64354 50.746, 4.46355 50.7551, 4.30794 50.7004, 3.93145 50.6895, 3.8798 50.7515, 3.68874 50.7747, 3.56821 50.7281, 3.17839 50.7561, 3.24515 50.713, 3.28845 50.5259, 3.47439 50.5337, 3.66415 50.4524, 3.65726 50.3686, 3.91011 50.3289, 4.02795 50.3584, 4.12439 50.2729, 4.23058 50.0693, 4.19482 49.9568, 4.44551 49.9372, 4.68472 49.9968, 4.70207 50.0956, 4.85042 50.1, 4.81917 49.9952, 4.89013 49.9089, 4.86309 49.8018, 5.26851 49.6966, 5.4112 49.6083, 5.4837 49.5075, 5.84181 49.5531, 5.90772 49.6391, 5.7357 49.8969, 5.89532 50.1122, 6.14766 50.1605, 6.17576 50.2354, 6.39938 50.3451, 6.3775 50.4396, 6.19691 50.5303, 6.19097 50.6397, 6.11071 50.7235, 5.89217 50.7552, 5.68206 50.7575))");

  CT ct;
  Polygon_with_holes_2 P;
  while(CGAL::IO::read_polygon_WKT(ifs, P)){
    const Polygon_2& poly = P.outer_boundary();
    ct.insert_constraint(poly);
    for(Polygon_with_holes_2::Hole_const_iterator it = P.holes_begin(); it != P.holes_end(); ++it){
      const Polygon_2& hole = *it;
      ct.insert_constraint(hole);
    }
  }

  for(Constraint_iterator cit = ct.constraints_begin();
      cit != ct.constraints_end();
      ++cit) {
    std::cout << "simplified polyline" << std::endl;
    for(Points_in_constraint_iterator vit =
          ct.points_in_constraint_begin(*cit);
        vit != ct.points_in_constraint_end(*cit);
        ++vit)
      std::cout << *vit << std::endl;
  }
  return 0;
}

I'm also compiling the code using CMake, with the following CMakeLists.txt file:

cmake_minimum_required(VERSION 3.1)

if (NOT CMAKE_BUILD_TYPE)
  set(CMAKE_BUILD_TYPE "Release" CACHE STRING "" FORCE)
endif()

project(demo)
find_package(CGAL)
add_executable(demo demo.cpp)
target_link_libraries(demo CGAL::CGAL)

This is an output that I got. I dumped the standard output of 2 runs into 2 text files and compared them using diff:

$ ./demo > out1.txt
$ ./demo > out2.txt
$ diff out1.txt out2.txt          
108a109,152
> 5.68787 50.8119
> 5.64166 50.8642
> 5.75882 50.9517
> 5.75839 51.035
> 5.85578 51.1446
> 5.24219 51.3052
> 5.07117 51.3935
> 4.66969 51.4264
> 4.53818 51.4824
> 4.38484 51.4493
> 4.43169 51.375
> 4.24205 51.354
> 4.07063 51.2507
> 3.88635 51.2002
> 3.59082 51.3044
> 3.44824 51.2415
> 3.36613 51.3693
> 3.07833 51.3027
> 2.5455 51.089
> 2.71077 50.8136
> 2.86315 50.7083
> 3.019 50.7736
> 3.17839 50.7561
> 3.56821 50.7281
> 3.68874 50.7747
> 3.8798 50.7515
> 3.93145 50.6895
> 4.30794 50.7004
> 4.46355 50.7551
> 4.64354 50.746
> 4.74861 50.8074
> 4.98168 50.7698
> 5.09584 50.7034
> 5.38958 50.748
> 5.47863 50.7235
> 5.68787 50.8119
> simplified polyline
> 4.24383 50.8196
> 4.29449 50.8892
> 4.43286 50.8947
> 4.48228 50.793
> 4.38314 50.7637
> 4.24383 50.8196
> simplified polyline
157,200d200
< simplified polyline
< 5.68787 50.8119
< 5.64166 50.8642
< 5.75882 50.9517
< 5.75839 51.035
< 5.85578 51.1446
< 5.24219 51.3052
< 5.07117 51.3935
< 4.66969 51.4264
< 4.53818 51.4824
< 4.38484 51.4493
< 4.43169 51.375
< 4.24205 51.354
< 4.07063 51.2507
< 3.88635 51.2002
< 3.59082 51.3044
< 3.44824 51.2415
< 3.36613 51.3693
< 3.07833 51.3027
< 2.5455 51.089
< 2.71077 50.8136
< 2.86315 50.7083
< 3.019 50.7736
< 3.17839 50.7561
< 3.56821 50.7281
< 3.68874 50.7747
< 3.8798 50.7515
< 3.93145 50.6895
< 4.30794 50.7004
< 4.46355 50.7551
< 4.64354 50.746
< 4.74861 50.8074
< 4.98168 50.7698
< 5.09584 50.7034
< 5.38958 50.748
< 5.47863 50.7235
< 5.68787 50.8119
< simplified polyline
< 4.24383 50.8196
< 4.29449 50.8892
< 4.43286 50.8947
< 4.48228 50.793
< 4.38314 50.7637
< 4.24383 50.8196

I tried perusing CGAL's documentation, but I couldn't find anything on this. Setting -DCMAKE_BUILD_TYPE=debug doesn't seem to make the iteration deterministic either. If it helps, I'm currently building the code on a M2 Mac running macOS 12.6, with the compiler AppleClang 14.0.0.14000029 and CGAL 5.5.1.

Answer 1

we have a pull request on github for making things in the Constrained_triangulation_plus_2 deterministic. github.com/CGAL/cgal/pull/8273

Answer 2

The documentation for Constrained_triangulation_plus_2 says that

The order of visit is undefined

You can make the order deterministic by sorting the points first. Something like this might work:

  set::set points;
  for(Constraint_iterator cit = ct.constraints_begin();
      cit != ct.constraints_end();
      ++cit) {
    std::cout << "simplified polyline" << std::endl;
    for(Points_in_constraint_iterator vit =
          ct.points_in_constraint_begin(*cit);
        vit != ct.points_in_constraint_end(*cit);
        ++vit)
      points.insert(*vit);
  }
  for(auto point const& : points) {
    std::cout << point << std::endl;
  }