Infinite loop when building intrinsic Delaunay triangulation in Heat_method_3

[oboes]

Issue Details

Using the Heat_method_3 package in Intrinsic_Delaunay mode with a mesh containing degenerate triangles will cause the algorithm to get stuck in an infinite loop.

The responsible loop can be found in internal/Intrinsic_Delaunay_triangulation_3.h:293. When flipping edges to build the Delaunay triangulation, the algorithm will try to add/remove the same cycle of edges to the stack indefinitely. (On the line just above there is also a int a = 0; declaration which serves no purpose).

This infinite loop is in turn caused by a division by zero when computing cotangents in internal/Intrinsic_Delaunay_triangulation_3.h:226. When encountering a degenerate triangle of the "needle" kind (having two vertices with identical positions), a NaN value will be produced, causing the is_edge_locally_delaunay function to return false, and finally the loop_over_edges to goes on forever.

(If the algorithm instead encounters a degenerate triangle of the "cap" kind, with 3 collinear but distinct points, then an inf value will be produced, the loop will terminate, and you will get a Eigen Decomposition in cotan failed assertion exception later on in the linear solver.)

The Heat_method_3 package make some common assumptions on the input mesh (such as it being a triangle mesh, etc) so it might be expected that it will fail if you feed it a mesh containing degeneracies, but as there is no assertions checking for that the software will just silently freeze at some point. I believe that degeneracies might also be created during the edge flipping phase due to floating-point rounding errors (this package only works with double data type so we can't use exact arithmetic there).

See the original paper describing the heat method for additional details on this algorithm.

Source Code

Minimal example :

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh/Surface_mesh.h>
#include <CGAL/Heat_method_3/Surface_mesh_geodesic_distances_3.h>

int main(int argc, char* argv[]) {

    using namespace CGAL;
    using Kernel = Exact_predicates_inexact_constructions_kernel;
    using Mesh = Surface_mesh<Kernel::Point_3>;

    Mesh mesh;
    Mesh::Vertex_index v1 = mesh.add_vertex( Mesh::Point(0, 0, 0) );
    Mesh::Vertex_index v2 = mesh.add_vertex( Mesh::Point(0, 0, 0) );
    Mesh::Vertex_index v3 = mesh.add_vertex( Mesh::Point(1, 0, 0) );
    Mesh::Vertex_index v4 = mesh.add_vertex( Mesh::Point(0, 1, 0) );
    mesh.add_face(v1, v2, v4);
    mesh.add_face(v2, v3, v4);
    mesh.add_face(v1, v3, v2);

    Mesh::Property_map<Mesh::Vertex_index, double> dmap
        = mesh.add_property_map<Mesh::Vertex_index, double>("v:dist").first;

    using Mode = Heat_method_3::Intrinsic_Delaunay;
    Heat_method_3::Surface_mesh_geodesic_distances_3<Mesh, Mode> heat(mesh);
    heat.add_source(v1);
    heat.estimate_geodesic_distances(dmap);

    return EXIT_SUCCESS; // this will never happen
}

Suggestions

I didn't submit a pull request for this issue because I am not sure of the preferred way to fix it. Below are some suggestions.

• Simply add an assertion checking for degeneracies and/or divisions by zero, so at least the algorithm will terminate.
• Catch NaN values and replace them by something sane which will produce the expected behaviour from the algorithm (if that's possible).
• Detect degenerate triangles and remove them from the intrinsic Delaunay triangulation, for example by using the remove_degenerate_faces function (but see also this related issue).
• Introduce some small random pertubation to remove degeneracies (after all the heat method is already a non-exact method for computing distances).
• Improve the Heat_method_3 package by implementing the method described in this recent paper by the same author which I believe could solve this and other issues related to non-manifold meshes (but that would require much more work).

Environment

• Operating system: Linux 5.8.10 x86_64
• Compiler: gcc 10.2.0
• Release or debug mode: both
• CGAL version: master/5.1
• Boost version: 1.72.0
• Eigen version: 3.3.7

[janetournois]

Hello @oboes for your detailed message.

I have made some experiments to try dealing with degenerate faces, and avoid manipulating NANs inside the code. I managed to avoid them but the infinite loop is still there.

For now I think that we should add a precondition that the input mesh should not have degenerate faces, and maybe refer to the mesh repair functions available in PMP. Maybe later we can think of an extension of the intrinsic triangulation that deals with degenerate faces, or implement the recent paper you mentioned.

@sloriot @afabri do you agree or have suggestions?

[sloriot]

+1

[oboes]

+1

[afabri]

I also fully agree. I had missed that the degenerate triangle was in the input, and had thought that the intrinsification had produced one. We still could ask Keenan Crane if he sees a trivial fix.

[janetournois]

Good. I'm preparing a PR for that and contacting Keenan Crane.

[oboes]

Today I read the mentioned paper N. Sharp & K. Crane 2020 and in fact the bulk of the article is about handling non-manifold input meshes by building some kind of double cover of the mesh (the "tufted cover") and then computing the Laplacian on it. Since the Surface_mesh class implements a half-edge data structure which anyway can't handle non-manifold edges (i.e. edges with at least 3 incident faces), we don't really need this tufted cover construction here.

If we only want to handle degeneracies then it is enough to implement the "mollification step" described in section 4.5 of the article, which is much less work, so I submitted the alternative PR #5037 to solve this issue. However it might require more tests to be sure that in meshes with no degeneracies, the end result is really the same with and without the mollification step.