CDT_plus_2 should prevents the insertion of twice the same constraint

[yobeonline]

Updated Issue Details, @lrineau 2018/08/02

The following test case should not loop, because it uses CDT_plus_2, that should prevent the insertion of twice the same constraint. See https://github.com/CGAL/cgal/issues/3267#issuecomment-409906745.

---

The old title was: [Probable infinite loop] Constrained 2D Delaunay triangulation : inserting a set of constraints takes forever

Issue Details

CGAL 4.11.3 introcuded a problem in the following code snippet (it was not present with previous releases (tested 4.11, 4.11.1, 4.11.2)). Inserting a constraint causes the insert_constraint method to never return.

The set of constraints are shown in the picture below:

Source Code

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/Constrained_triangulation_plus_2.h>
#include <CGAL/Triangulation_face_base_2.h>
#include <CGAL/Triangulation_vertex_base_2.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/Projection_traits_xy_3.h>

using K1 = CGAL::Exact_predicates_inexact_constructions_kernel;
using K = CGAL::Projection_traits_xy_3<K1>;

using Vb2 = CGAL::Triangulation_vertex_base_2<K>;
using Fb2 = CGAL::Constrained_triangulation_face_base_2<K>;

using TDS2 = CGAL::Triangulation_data_structure_2<Vb2,Fb2>;
using Itag2 = CGAL::Exact_predicates_tag;
using CDT2_Base = CGAL::Constrained_Delaunay_triangulation_2<K,TDS2,Itag2>;

using Triangulation = CGAL::Constrained_triangulation_plus_2<CDT2_Base>;
using point = Triangulation::Point;

void main(void)
{
  Triangulation t;

  point const p1(866863.099976, 2025662.79999,0);
  point const p2(866876.599976, 2025663.20001,0);
  point const p3(866873.799988, 2025719,0);
  point const p4(866859.599976, 2025718.09998,0);

  point const p6(866872.099976, 2025645.29999,0);
  point const p7(866888.599976, 2025646.20001,0);
  point const p8(866887.599976, 2025663.59998,0);
  point const p9(866870.599976, 2025663,0);

  auto const insert = [&t,p1,p2,p3,p4,p6,p7,p8,p9]()
  {
    { // green
      t.insert_constraint(p1,p2);
      t.insert_constraint(p2,p3);
      t.insert_constraint(p3,p4);
      t.insert_constraint(p4,p1);
    }

    { // blue
      t.insert_constraint(p6,p7);
      t.insert_constraint(p7,p8);
      t.insert_constraint(p8,p9); // *
      t.insert_constraint(p9,p6);
    }
  };

  insert();
  insert(); // this call never ends beacuse of the marked line in the insert lambda.
}  

I admit the use case is a bit contrived (inserting twice the same data). I don't know how to generate it in a simpler manner. I figured it was still worth submitting it.

Environment

• Windows 64 bits
• Compiler: vs 2013 32 and 64
• Release or debug mode: both
• Specific flags used (if any): none
• CGAL version: 4.11.3
• Boost version: 1.65.1
• Other libraries versions if used (Eigen, TBB, etc.): none

[lrineau]

Thanks for the isolation of that issue with a small subset of data.

That issue here again is the cascading of construction, using an inexact arithmetic. Here by inserting the segments several time, you just worsen the situation.

The good news is that the patch introduced in CGAL-4.11.2 and CGAL-4.11.3 can be reused. By default, it snap newly constructed points to existing points with a distance of 4 ULP (ignoring between 2 and 3 bits of precision). You can add the following line:

#define CGAL_CDT_2_INTERSECTION_SNAPPING_ULP_DISTANCE 8

at the top of your .cpp file, and it will use that distance 8 instead of 4, ignoring between 3 and 4 bits of precision. And that fixes your test case. Maybe with your real data you will need a bigger value. Each time you double that value, you ignore one more bit of precision (among the 53 bits of precision of the double type).

[yobeonline]

Thank you for the comment, this define may very well make my day!

(I don't known where I should have the following discussion)

The data I have is rounded to 10e^-3 (although it is not the case here) because the rest of the application will not support faces that are too small (area-wise). So maybe I should use this define to tell cgal that I don't mind it snapping far away (with regard to 4 or even 8 ulps).

I was browsing the documentation to see if triangulating with a fixed point type would be the more natural solution but I am not sure I understood the kernel I should use : Cartesian<int>?

[lrineau]

@yobeonline commented on Aug 1, 2018, 11:01 AM GMT+2:

The data I have is rounded to 10e^-3 (although it is not the case here) because the rest of the application will not support faces that are too small (area-wise). So maybe I should use this define to tell cgal that I don't mind it snapping far away (with regard to 4 or even 8 ulps).

Be careful that 10^-3 does not give a relative precision. If your range of data coordinates is between 0. and 1., that is 3 digits (in base 10) of precision. But if your data coordinates are between 100000 and 999999 then that mean 9 digits of precision. You can use log2 to convert to bits of precision (in base 2). For example 9 digits of precision is about 30 bits of precision. As double has 53 bits of precision, you can waste 23 of them. And that mean, if I am not wrong, that you can use:

#define CGAL_CDT_2_INTERSECTION_SNAPPING_ULP_DISTANCE 4194304 // 2^22

I have never tried the code with such a big value.

I was browsing the documentation to see if triangulating with a fixed point type would be the more natural solution but I am not sure I understood the kernel I should use : Cartesian ?

If I were you, I would not try that approach. The kernel CGAL::Exact_predicates__inexact_constructions_kernel is very different from Cartesian<double> or Simple_cartesian<double> (the later is without reference counting). The difference is that its predicates, such as the Orientation_2 predicates, has a layer of filters, and exact arithmetics, to be both fast (in most cases) and exact (with degenerate cases). If you start with something like Simple_cartesian<int>, you will have two kind of issues:

1. First, note that int, on Windows with Visual Studio, is a 32 bits integer (both in 32 and 64 bits compilations). That means that you will have 32 bits of precision. That seems enough, but unless floating points that can scale from -2^-1022 to a bit more than 2^1023, you also need extra bits when you multiply two integers. Predicates in 2D triangulation are of the second degree, as far as I know, so that is bearable. I would anyway choose long long int instead of int, to have 64 bits in both 32 bits and 64 bits compilations.

2. The difficulty will be to get the layers of filters and exactness of the CGAL::Exact_predicates__inexact_constructions_kernel kernel. The exact layer would also require that you can convert long long int to the exact type (CGAL::Gmpq), and that is easy. The Interval_nt filtering would also probably work. But I do not know how to adapt the static filters, that are specific to the number type double. You could drop the static filters but that would impact performances.

My advice is to keep the CGAL::Exact_predicates__inexact_constructions_kernel kernel.

[yobeonline]

1. I understood the "relative issue" of the 10e-3 rounding. I will give it a shot a report the results here.
2. Thank you for the detailed explanation of kernels. I was not particularly hyped by the idea of trying another one, so if you think I should not, well, so be it!

[yobeonline]

Hi again, I was wondering why did you close this ticket? I agree that changing the snapping distance resolves this occurence of the problem, but I am concerned with the fact that a call to insert_constraint may never return. Unless you can confirm that my example will eventually return, I suggest you reopen this ticket.

[lrineau]

I'll write a big comment soon.

For the moment, I want to say that I was able to reduce the test case that loops to:

#define CGAL_CDT_2_INTERSECTION_SNAPPING_ULP_DISTANCE 7 // 4194304
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/Constrained_triangulation_plus_2.h>
#include <CGAL/Triangulation_face_base_2.h>
#include <CGAL/Triangulation_vertex_base_2.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/Projection_traits_xy_3.h>

using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using Fb2 = CGAL::Constrained_triangulation_face_base_2<K>;
using Vb2 = CGAL::Triangulation_vertex_base_2<K>;

using TDS2 = CGAL::Triangulation_data_structure_2<Vb2,Fb2>;
using Itag2 = CGAL::Exact_predicates_tag;
using CDT2_Base = CGAL::Constrained_Delaunay_triangulation_2<K,TDS2,Itag2>;

using Triangulation = CGAL::Constrained_triangulation_plus_2<CDT2_Base>;
using point = Triangulation::Point;

int main()
{
  Triangulation t;

  point const p1(866863.099976, 2025662.79999);
  point const p2(866876.599976, 2025663.20001);
  point const p3(866873.799988, 2025719);

  point const p8(866887.599976, 2025663.59998);
  point const p9(866870.599976, 2025663);

  t.insert_constraint(p1,p2);
  t.insert_constraint(p2,p3);
  t.insert_constraint(p8,p9);
  t.insert_constraint(p1,p2);
  t.insert_constraint(p8,p9);
}

See also https://gist.github.com/2c4882d4fc67f3c88bd9bfa0141093b2

[lrineau]

I can reduce the issue to:

int main()
{
  Triangulation t;

  point const p1( 866874.52509081131, 2025663.1385288462 ); // p26
  point const p2( 866874.52509079711, 2025663.1385288464 ); // p40
  point const p3( 866874.52509081364, 2025663.1385288467 ); // p45
  point const p4( 866874.52509080991, 2025663.1385288469 ); // p50
  point const p5( 866874.52509081631, 2025663.1385288464 ); // p52

  t.insert_constraint(p1,p3);
  t.insert_constraint(p3,p2);
  t.insert_constraint(p4,p5);
}

The point number in the comments are the numbering of points created by the test case of my last comment.

With this patch to CGAL, I can compile and run the following modified program:

#define CGAL_CDT_2_DEBUG_INTERSECTIONS 1
#define CGAL_CDT_2_INTERSECTION_SNAPPING_ULP_DISTANCE 7 // 4194304
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/Constrained_triangulation_plus_2.h>
#include <CGAL/Triangulation_face_base_2.h>
#include <CGAL/Triangulation_vertex_base_2.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/Projection_traits_xy_3.h>

using K1 = CGAL::Exact_predicates_inexact_constructions_kernel;
using K = K1; // CGAL::Projection_traits_xy_3<K1>;

template <typename Vb>
class My_vertex_base : public Vb {
  std::size_t time_stamp_;
public:
  My_vertex_base() : Vb(), time_stamp_(-1) {
    // std::cerr << "My_vertex_base(" << (void*)(this) << ") time_stamp_ is "
    //           << this->time_stamp_ << "\n";
  }

  My_vertex_base(const My_vertex_base& other) :
    Vb(other),
    time_stamp_(other.time_stamp_)
  {}

  typedef CGAL::Tag_true Has_timestamp;

  std::size_t time_stamp() const {
    return time_stamp_;
  }
  void set_time_stamp(const std::size_t& ts) {
    time_stamp_ = ts;
  }

  template < class TDS >
  struct Rebind_TDS {
    typedef typename Vb::template Rebind_TDS<TDS>::Other Vb2;
    typedef My_vertex_base<Vb2> Other;
  };
};

#ifdef CGAL_CDT_2_DEBUG_INTERSECTIONS
using Vb2 = My_vertex_base<CGAL::Triangulation_vertex_base_2<K>>;
#else
using Vb2 = CGAL::Triangulation_vertex_base_2<K>;
#endif

using Fb2 = CGAL::Constrained_triangulation_face_base_2<K>;

using TDS2 = CGAL::Triangulation_data_structure_2<Vb2,Fb2>;
using Itag2 = CGAL::Exact_predicates_tag;
using CDT2_Base = CGAL::Constrained_Delaunay_triangulation_2<K,TDS2,Itag2>;

using Triangulation = CGAL::Constrained_triangulation_plus_2<CDT2_Base>;
using point = Triangulation::Point;

int main()
{
  std::cerr.precision(17);
  Triangulation t;

  point const p1( 866874.52509081131, 2025663.1385288462 ); // p26
  point const p2( 866874.52509079711, 2025663.1385288464 ); // p40
  point const p3( 866874.52509081364, 2025663.1385288467 ); // p45
  point const p4( 866874.52509080991, 2025663.1385288469 ); // p50
  point const p5( 866874.52509081631, 2025663.1385288464 ); // p52

  t.insert_constraint(p1,p3);
  t.insert_constraint(p3,p2);
  t.insert_constraint(p4,p5);
}

and its execution gives:

CDT_p_2::   insert_constraint( #1= 866874.52509081131 2025663.1385288462 , #2= 866874.52509081364 2025663.1385288467 )
CDT_p_2::insert_subconstraint( #1= 866874.52509081131 2025663.1385288462 , #2= 866874.52509081364 2025663.1385288467 )
CDT_p_2::insert_subconstraint stack( #1= 866874.52509081131 2025663.1385288462 , #2= 866874.52509081364 2025663.1385288467 )
CDT_p_2::   insert_constraint( #2= 866874.52509081364 2025663.1385288467 , #3= 866874.52509079711 2025663.1385288464 )
CDT_p_2::insert_subconstraint( #2= 866874.52509081364 2025663.1385288467 , #3= 866874.52509079711 2025663.1385288464 )
CDT_p_2::insert_subconstraint stack( #2= 866874.52509081364 2025663.1385288467 , #3= 866874.52509079711 2025663.1385288464 )
CDT_p_2::   insert_constraint( #4= 866874.52509080991 2025663.1385288469 , #5= 866874.52509081631 2025663.1385288464 )
CDT_p_2::insert_subconstraint( #4= 866874.52509080991 2025663.1385288469 , #5= 866874.52509081631 2025663.1385288464 )
CDT_p_2::insert_subconstraint stack( #4= 866874.52509080991 2025663.1385288469 , #5= 866874.52509081631 2025663.1385288464 )
  CT_2::intersection: intersection SNAPPED to an existing point 866874.52509081631 2025663.1385288464
  intersection (#4, #5) x (#2, #3) true  vertex #5 = 866874.52509081631 2025663.1385288464
P_c_h::split-constraint( #2= 866874.52509081364 2025663.1385288467 , #3= 866874.52509079711 2025663.1385288464 , #5= 866874.52509081631 2025663.1385288464 )
CDT_p_2::insert_subconstraint( #2= 866874.52509081364 2025663.1385288467 , #5= 866874.52509081631 2025663.1385288464 )
CDT_p_2::insert_subconstraint stack( #2= 866874.52509081364 2025663.1385288467 , #5= 866874.52509081631 2025663.1385288464 )
CDT_p_2::insert_subconstraint( #5= 866874.52509081631 2025663.1385288464 , #3= 866874.52509079711 2025663.1385288464 )
CDT_p_2::insert_subconstraint stack( #5= 866874.52509081631 2025663.1385288464 , #3= 866874.52509079711 2025663.1385288464 )
  CT_2::intersection: intersection SNAPPED to an existing point 866874.52509081364 2025663.1385288467
  intersection (#5, #3) x (#1, #2) true  vertex #2 = 866874.52509081364 2025663.1385288467
CDT_p_2::insert_subconstraint( #1= 866874.52509081131 2025663.1385288462 , #2= 866874.52509081364 2025663.1385288467 )
CDT_p_2::insert_subconstraint stack( #1= 866874.52509081131 2025663.1385288462 , #2= 866874.52509081364 2025663.1385288467 )
  intersection at #2= 866874.52509081364 2025663.1385288467
P_c_h::split-constraint( #5= 866874.52509081631 2025663.1385288464 , #3= 866874.52509079711 2025663.1385288464 , #2= 866874.52509081364 2025663.1385288467 )
CDT_p_2::insert_subconstraint stack( #2= 866874.52509081364 2025663.1385288467 , #3= 866874.52509079711 2025663.1385288464 )
CDT_p_2::insert_subconstraint stack( #5= 866874.52509081631 2025663.1385288464 , #2= 866874.52509081364 2025663.1385288467 )
  intersection at #5= 866874.52509081631 2025663.1385288464
CDT_p_2::insert_subconstraint stack( #4= 866874.52509080991 2025663.1385288469 , #5= 866874.52509081631 2025663.1385288464 )
  CT_2::intersection: intersection SNAPPED to an existing point 866874.52509081631 2025663.1385288464
  intersection (#4, #5) x (#2, #3) true  vertex #5 = 866874.52509081631 2025663.1385288464
P_c_h::split-constraint( #2= 866874.52509081364 2025663.1385288467 , #3= 866874.52509079711 2025663.1385288464 , #5= 866874.52509081631 2025663.1385288464 )

then it loops.

[lrineau]

I have closed the ticket, because there is not much we can do with input segments that intersect, and an inexact construction of intersection, using floating-point numbers. We can find workarounds, like that trick with the snapping up to a distance, but, as you can see with the example with three segments in the previous comment, that workaround has limits: all the intersections were snapped, and that is why that looped.

There is one thing that we could do, in your specific test case that insert constraints twice, is to prevent the insertion of twice the same constraint.

@afabri, there is that piece of code in CGAL/Constrained_triangulation_plus_2.h: https://github.com/CGAL/cgal/blob/815e23538020280f5102177181a5e7297af6e88f/Triangulation_2/include/CGAL/Constrained_triangulation_plus_2.h#L267-L278 but the hierarchy.insert_constraint(va, vb) never returns NULL: https://github.com/CGAL/cgal/blob/815e23538020280f5102177181a5e7297af6e88f/Triangulation_2/include/CGAL/Polyline_constraint_hierarchy_2.h#L865-L891

[lrineau]

I reopen the issue. I have modified its title, and amended the first message as well.

[afabri]

The problem is certainly more general, than identical segments. The same probably happens when you have overlapping collinear segments.

[lrineau]

I agree, but either you remove the test, or you fix it.

[lrineau]

@lrineau commented on Aug 2, 2018, 2:19 PM GMT+2:

There is one thing that we could do, in your specific test case that insert constraints twice, is to prevent the insertion of twice the same constraint.

I have implemented that in the new PR #3271. But there will not be a bug-fix release for the CGAL-4.11-branch any more. The patch will be in CGAL-4.12.1, probably late August once I come back from vacation.

[lrineau]

Fixed by #3271.