CGAL Polyline Simplification results in self-intersection

Question

I'm currently having a bit of trouble with CGAL's Polyline Simplification.

More specifically, for the following example, PS::simplify(ct, Cost(), Stop(0.2)) results in a self-intersecting polyline. In the image below, the blue polyline is the input polyline into PS::simplify() while the green polyline is the resulting (output) polyline. The red arrow points to the self-intersection in the resulting polyline.

Further below, I have copied and pasted my code from 2 files simplify_test.cpp and CMakeLists.txt. With the required libraries installed, to run this example, you may place them in the same directory, cd to that directory, and run the following in your terminal:

$ cmake .
$ make
$ ./simplify_test

This outputs the resulting coordinates (green polyline in the image below), which contain the self-intersection.

I would like to know:
(1) why this results in a self-intersection
and (2) what can be done so that it does not result in a self-intersection.

Thank you for your help!

Here is the simplification code in a file named simplify_test.cpp:

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polyline_simplification_2/simplify.h>

#include <iostream>
#include <string>
#include <vector>

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Point_2<K> Point;
typedef std::vector<Point> Polyline;
namespace PS = CGAL::Polyline_simplification_2;
typedef PS::Vertex_base_2<K> Vb;
typedef CGAL::Constrained_triangulation_face_base_2<K> Fb;
typedef CGAL::Triangulation_data_structure_2<Vb, Fb> TDS;
typedef CGAL::Constrained_Delaunay_triangulation_2<K, TDS, CGAL::Exact_predicates_tag> CDT;
typedef CGAL::Constrained_triangulation_plus_2<CDT> CT;
typedef PS::Stop_below_count_ratio_threshold Stop;
typedef PS::Squared_distance_cost Cost;

void print_coords(Polyline polyline) {
  std::cout << std::endl;
  std::cout << "Simplified coordinates:" << std::endl << std::endl;

  // Print out the simplified coordinates
  unsigned int i = 0;
  for (Point coord : polyline) {
    std::cout << "[ ";
    std::cout << coord.x();
    std::cout << ", ";
    std::cout << coord.y();
    std::cout << " ]";
    if (i != polyline.size() - 1) std::cout << ",";
    i++;
  }
  std::cout << std::endl << std::endl;
}

void simplify_test() {
  // Hard code a minimum working example where running PS::simplify results in
  // self-intersections. There are no self-intersections when {27, 9} is
  // omitted.
  std::vector<std::vector<int>> coords = {
      {64, 20}, {33, 27}, {27, 9}, {33, 18}, {44, 18}, {44, 8},
      {24, 0},  {0, 13},  {9, 49}, {84, 41}, {83, 29}, {64, 20},
  };
  // Simplification outputs:
  // [ 64, 20 ],[ 27, 9 ],[ 44, 18 ],[ 24, 0 ],
  // [ 9, 49 ],[ 83, 29 ],[ 64, 20 ],[ 64, 20 ]

  // Create polyline for simplifying later
  Polyline polyline;

  // Insert coordinates into polyline
  for (std::vector<int> coord : coords) {
    Point pt(coord[0], coord[1]);
    polyline.push_back(pt);
  }

  // Insert polyline into ct and run simplify()
  CT ct;
  ct.insert_constraint(polyline.begin(), polyline.end());
  PS::simplify(ct, Cost(), Stop(0.2));
  Polyline polyline_simplified;

  // Transfer simplified coordinates from ct to polyline for easy handling
  auto cit = ct.constraints_begin();
  for (auto vit = ct.points_in_constraint_begin(*cit);
       vit != ct.points_in_constraint_end(*cit); vit++) {
    polyline_simplified.push_back(*vit);
  }

  print_coords(polyline_simplified);
}

int main() {
  simplify_test();
}

Here is the CMakeLists.txt file.

cmake_minimum_required(VERSION 3.1)
project(simplify_test LANGUAGES CXX)

set(CMAKE_CXX_STANDARD 20)
set(CMAKE_BUILD_TYPE Release)

# Detecting appropriate compiler
if (APPLE)
  set(CMAKE_C_COMPILER "/usr/local/opt/llvm/bin/clang")
  set(CMAKE_CXX_COMPILER "/usr/local/opt/llvm/bin/clang++")
elseif(UNIX) # implicit AND NOT APPLE
  set(CMAKE_CXX_COMPILER "g++-10")
  set(CMAKE_C_COMPILER "gcc-10")
endif()

# Finding appropriate packages
find_package(CGAL)

# Adding executables needed
add_executable(simplify_test ./simplify_test.cpp)

# Linking appropriate libraries required
target_link_libraries(
  simplify_test
  CGAL::CGAL
)

Answer 1

In the C++ code below, I essentially replaced Polyline with Polygon_with_holes_2. For all values of stop, I now get topologically valid simplified polygons. I edited your output function. You can copy and paste the output from the new function print_coords_for_geogebra() directly into Geogebra.

Here is the edited version of simplify_test.cpp. You can compile it with the same CMakeLists.txt that you included in your original post.

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polyline_simplification_2/simplify.h>
#include <CGAL/Polygon_with_holes_2.h>
#include <iostream>
#include <vector>

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Point_2<K> Point;
namespace PS = CGAL::Polyline_simplification_2;
typedef PS::Stop_below_count_ratio_threshold Stop;
typedef PS::Squared_distance_cost Cost;
typedef CGAL::Polygon_2<K> Polygon;
typedef CGAL::Polygon_with_holes_2<K> Polygon_with_holes;

void print_coords_for_geogebra(Polygon ext_ring)
{
  unsigned int i = 0;
  std::cout << "Polygon(";
  for (Point coord : ext_ring) {
    std::cout << "("
              << coord.x()
              << ", "
              << coord.y()
              << ")";
    if (i != ext_ring.size() - 1) std::cout << ", ";
    i++;
  }
  std::cout << ")" << std::endl << std::endl;
}

void simplify_test2(Polygon_with_holes polygon, double stop)
{
  Cost cost;
  polygon = PS::simplify(polygon, cost, Stop(stop));
  print_coords_for_geogebra(polygon.outer_boundary());
}

int main(int argc, char* argv[])
{
  std::vector<std::vector<int> > coords = {
    {64, 20}, {33, 27}, {27, 9}, {33, 18}, {44, 18}, {44, 8},
    {24, 0},  {0, 13},  {9, 49}, {84, 41}, {83, 29}, {64, 20},
  };

  // Insert coordinates into the external ring of a polygon with holes
  Polygon ext_ring;
  for (std::vector<int> coord : coords) {
    Point pt(coord[0], coord[1]);
    ext_ring.push_back(pt);
  }
  Polygon_with_holes polygon(ext_ring);
  print_coords_for_geogebra(polygon.outer_boundary());
  for (double stop = 1.0; stop > 0.0; stop -= 0.1) {
    simplify_test2(polygon, stop);
  }
  return EXIT_SUCCESS;
}

Answer 2

I further reduced it to a polyline std::vector<std::vector<int> > coords = { {64, 20}, {33, 27}, {27, 9}, {33, 18}, {44, 18}, {24, 0} };

You found a bug in the points in constraint iterator. As a workaround you should use the Vertex_in_constraint_iterator and then call the method point(). I will file an issue on github and a fix.