A failure case in constrained Delaunay triangulation in mesh intersection

[BrunoLevy]

There is an example here for which the Constrained Delaunay Triangulation generates two different triangulations for the same convex quadrilateral.

analysis Interval filtering deactivated FPE activated (but was not triggered)

• with CDT2dBase::exact_incircle_ = true and exact_nt assertion fail in orient2d_lifted_SOS(), Delta3_sign is ZERO
• with CDT2dBase::exact_incircle_ = false and exact_nt, infinite loop in Delaunayize_new_edges()
• with CDT2dBase::exact_incircle_ = true and expansion_nt, no error detected, but wrong result
• with CDT2dBase::exact_incircle_ = false and expansion_nt, no error detected, but wrong result
  • debug mode, with exact_incircle_=true and exact_nt, assertion fail edge_is_delaunay

I get a different result with this dataset with expansion_nt (correct result) and exact_nt (wrong result). For both, using exact_incircle_ = false and no filters.

With exact_incircle and PR22, during the first triangle, we have an assertion fail (edge is not Delaunay) right after inserting the first constraint. The first triangle has three vertices to insert, they are aligned, and there are three constraints connecting them (one of the constraints is superfluous !). Aha, Delaunay condition is not satisfied before inserting the constraint. The culprit seems to be insertion of a point in an edge. So we insert a vertex on a border edge then Delaunay condition is not satisfied, how can this happen ? Now I suspect that the two triangles that I enqueue in the queue of triangles to be Delaunayized are not rotated in such a way that the newly inserted vertex is vertex 0, but I am unsure, because I also observed the bug with the naive implementation that does not have this requirement, but is this the case ? Let us check with the naive implementation -> has the same problem. So it may be my incircle that is wrong, but I doubt it. Or maybe it is my consistency check that is wrong ? Little spoon debugging Let us see which pair of triangles violate the Delaunay condition. Well, in fact, of course I cannot flip the two triangles generated by inserting a point on an edge, so I need to re-read the definition to see how my consistency check should be implemented (seems that my consistency check is bugged)

Ooo, I obtain a different result when using expansion_nt and exact_nt in orient_2dlifted_SOS_projected(): on three_cubes_4_faces.obj, correct result with expansion_nt, assertion fail Tedge_is_Delaunay() with exact_nt. But with three_cubes.obj, both expansion_nt and exact_nt fail.

hypothesis So there is a difference between expansion_nt and exact_nt. It is probably expansion_nt that is wrong. The Delta3_sign==ZERO in orient2d_lifted_SOS() gives me the impression that there are two duplicated points that were not identified as such. I'm also suspecting insertion of points on the border of the triangulation. In that case, are we able to get all the triangles to Delaunayize ?

debugging strategy

• [x] Implement stored coordinates as exact_nt to see whether it makes a difference.
• [x] Store an exact_nt in each expansion_nt and compare them after each operation to find the faulty operation. (not needed, the problem was different, see below)
• [x] Debug with the "little spoon", and generate outputs for all steps in my simple example with 6 triangles.
• [x] Display the points that trigger a ZERO determinant in the predicate

[BrunoLevy]

The two faulty triangles form a symmetric shape: in_circle() should return the same value for both configurations, the four points are cocyclic (but they are not detected as such by my predicate). Maybe it comes from the inexact handling of the "lifted coordinates". I can try computing them exactly in the exact_nt version. (OK, but this still does not explain the difference between expansion_nt and exact_nt, until there was an FPE, but I think I already checked...)

[BrunoLevy]

OK, so now I have an exact_nt class (that will not suffer from overflows/underflows). I had a bug in the compare() function (on the easy case with different signs, silly me !), now fixed. Now it seems I'm able to reproduce all the results of the expansion_nt-based pipeline (including the bug !, but it will be easier to analyze...)

[BrunoLevy]

• Implemented incircle() with exact evaluation of the lifting coordinates x^2+y^2 with exact_nt
• Sometimes we have infinite flipflop loops in Delaunayize_new_edges(). It means both in_circle(v1,v2,v3,v4) and in_circle(v1,v2,v4,v3) return POSITIVE (impossible, so the bug is probably there)
• OK, the in_circle() predicate is the culprit ! It says that both the flipped and non-flipped configurations are non-Delaunay, so it flips forever.

Found one of the problems: operator-(vec3Hex, vec3Hex) was reversed ! (always the same confusion between make_vector(p1,p2) = p2-p1 and operator-(p1,p2)

Now I mostly get the same results as with expansion_nt. In 7sins.obj, I get two quads not triangulated the same. Investigating that...

• OK, assertion fails edge_is_Delaunay()

[BrunoLevy]

It seems I got it (yeehaa !)

• one needs to take care when to Delaunayize vertex neighbors when handling intersecting constraints. If done after flipping the edges to enforce the constraint, then the connectedness of the conflict zone may be destroyed, preventing Delaunay condition to be recovered for all triangles
• the incircle2d(p1,p2,p3,p4) predicate in its usual form, that returns POSITIVE if p4 is in the circumscribed circle of triangle p1,p2,p3, supposes that p1,p2,p3 is oriented counterclockwise, which is not always the case in the intersection code (that projected the input triangle along its less-varying direction, without checking for orientation). The rest of the code in CDT2d was taking into account the global orientation for orient2d(), but not for incircle2d() !!
• incircle2d() with exact_nt: rethought computation and symbolic perturbation. Computation done with translated version of predicate (p4 -> origin), and symbolic perturbation using not-translated version (that has the original vertex lengths, and perturbations directly appear)
• incircle2d() with expansion_nt: rethought computation (symbolic perturbation remains the same): everything done with not-translated form, so that everything is in terms of original vertex lengths. The little game with w coordinates is simpler, because (approximated) lengths do not have a w coordinate.
• more care needs to be taken with approximated lengths: an additional quad convexity test using orient2d() is needed to avoid creating degenerate triangles or even slightly flipped triangles (that may also appear temporarily), this is the exact_incircle flag in CDT2d (set to false for expansions/approximated lengths)

[BrunoLevy]

done Post-debugging refactoring TODO list:

• [x] cleanup, simplifications in exact_nt predicates I'd like to be able to write cleaner code for SOS, like:

  Sign s = sos(
     p0, { return ... },
     p1, { return ...},
     p2, { return ...},
     p3, { return ...}
     } 

• [x] interval filters for exact_geometry implemented with exact_nt
• [x] incircle2d_SOS(), projected_incircle2d_SOS()
• [x] generic rationalg and vec2Hg and vec3Hg templates
• [x] test instanciations rationalg, rationalg, vec2Hg, vec2Hg
• [x ] test rational_nt::simplify()
• [x ] vec2Hex::simplify() and faster vec3HEx lexico comparison, canonicalized vec3HEx for point-to-vertex map

[BrunoLevy]

Still a problem with Spiral_sphere_ornament_2013-11-23.stl (from Thingy10K), assertion fail with inserted point outside boundary (!!!)

May come from a super skinny triangle nearly identical to a segment oriented towards (1,1,1) (worst case for chosing a projection direction).

-> filed new issue for this one (#110)

[BrunoLevy]

fixed Still a problem with PR21/bug_CDT2d.obj and PR21/bug_CDT2d2.obj (both assertion fail Delta_sign != ZERO in side3 predicate -> reactivated inexact incircle mode in CDT2d .

[BrunoLevy]

fixed A problem with DATA/Small/test_isect.obj that only appears under Valgrind:

valgrind bin/intersect test_isect.obj 

Result: Assertion failed: orient2d(v1,v2,v3) != ZERO (got also Assertion failed: nb_z != 3)

Fortunately, it also appears with sys:multithreading=false (will be easier to debug/reproduce), but the fact that it only happens under Valgrind is annoying...

It seems that valgrind does not support changing the rounding mode !

Problem does not appear with typedef intervalRN interval_nt (instead of intervalRU).

TODO: check difference of performance between intervalRN and intervalRU to see whether RU is worth the additional pain. If RU is significantly faster, we will need to dynamically detect Valgrind (e.g., by testing a rounding operation) and deactivate interval filtering if rounding is not correct.

[BrunoLevy]

Problems were:

• delaunayize_vertex_neighbors() needs to be called before anything else when creating a new vertex from constraints intersections
• in_sphere() predicate with (imprecise but coherent) precomputed lifted coordinates may create triangle flips, so an additional check for convex quadrilateral is needed in CDT2d
• symbolic perturbations were probably not correct (rewrote them in a much easier form)