In mesh_surface_intersection, point outside boundary assertion failure in constrained Delaunay triangulation

[BrunoLevy]

With Spiral_sphere_ornament_2013-11-23.stl (from Thingy10K), in debug mode:

  $ bin/intersect Spiral_sphere_ornament_2013-11-23.stl

assertion fail with inserted point outside boundary (!!!)

[BrunoLevy]

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

• Check that max normal component is non-zero: yes (pb is not there)
• Check whether the problem happens at the same time with exact_nt and expansion_nt: yes (pb is not there)
• Replace locate() with naive "scan all triangles" implementation: same story ! (pb is not in optimized locate)
• FP underflow / overflow in orient3d() ? yes ! (in interval arithmetics, filter for orient_2d and incircle_2d_SOS) Deactivate filter and see what happens: still get point outside boundary !!
• Let's see where this point comes from (which triangles, what do they look like ?)
  • we are in triangle 4079
  • the other triangle is 25590
  • the point to be inserted is on 4089.E2 and 25590.P2
• Could it be a bug in triangle_intersections() ?
• The 2D triangulation right before inserting the vertex that crashes contains 4 vertices, two triangles. When looking at it, it looks like a single segment, with three points at the same location. One of the triangles looks like a segment, the other one looks like a single point. Let's double check for zero-area triangles: there is no zero-area triangle
• Extracted the two triangles from Spiral_sphere_ornament_2013-11-23.stl that cause the problem, posted the file here (Experiment/DATA/PR23)

OK, so I understand nothing, let's go "little spoon debugging mode", we only got five points and two triangles, so let us double check the result of orient_3d() and the 2D cases. Another possibility is a bug in orient3d() (a bit unlikely but who knows). So we could re-implement it with exact_nt and see what happens...

The two faulty triangles are like that: one of them has two nearly colocated vertices, and is sharing a vertex with the other one (that has a rather normal shape). Even by zooming a lot on the two nearly colocated vertices, one cannot see that they are not the same, but in CDT2d, projected_orient2d() never returned ZERO, so the triangle is probably valid (but very very degenerate).

When displaying the combinatorial information in the triangle-triangle intersection, it seems that there are incoherencies.

Next step: do a little drawing, with the points and edges labeled, to figure out what the combinatorial information says.

Activated debug message in triangle_triangle_intersection()

Problem may be in triangle_intersection(), intersect_edge_edge_1d() that checks for vertex-in-edge in 1d. nope (was not that).

Gotcha! The function uses an inexact Geom::triangle_axis() function that may choose not the best projection axis. And the orient_2d() != ZERO checks I have done were not there ! Tested: with exact triangle axis computation it works.