3D Mesher with SDF Input does not terminate

stefango74 commented 3 months ago

Please use the following template to help us solving your issue.

Issue Details

CGAL 5.6. 3D Mesher used on an SDF will not terminate and remains in an infinite loop. This is not the specified behavior stated in the user manual (guaranteed to terminate for facet_angle = 30.0 degrees)

Source Code

Here is the code fragment that calls the mesher: Labeled_mesh_domain ldomain = Labeled_mesh_domain::create_implicit_mesh_domain(sdfInitial, Bbox_3(-1.0, -1.0, -1.0, 1.0, 1.0, 1.0)); // [-1,1] cube bounding the SDF Mesh_criteria_3 lcriteria(facet_angle = 30.0, edge_size = 0.1); LabeledC3t3 c3t3 = make_mesh_3(ldomain, lcriteria, manifold());

sdfInitial() returns the signed distance function on the [-1, 1] 3D grid, for points outside projected onto its bounding box

Environment

Operating system (Windows/Mac/Linux, 32/64 bits): Ubuntu 22.04
Compiler: g++ (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0
Release or debug mode: Debug mode
Specific flags used (if any): g++ -O0 -g3 -Wall -c -fmessage-length=0
CGAL version: 5.6
Boost version: 1.74
Other libraries versions if used (Eigen, TBB, etc.): g++ -lgmp -lz -ltbbmalloc -ltbb -lmpfr (tbb is disabled in compilation though)

stefango74 commented 3 months ago

I think I have located the reason for the infinite loop (but not the problem inside the mesher which causes it)

Specifically, the infinite loop where it does not exit anymore, is in Labeled_mesh_domain_3.h: // Else lets find a point (by bisection) // Bisection ends when the point is near than error bound from surface while(true) { // If the two points are enough close, then we return midpoint if ( squared_distance(p1, p2) < rdomain.squared_errorbound ) { CGAL_assertion(value_at_p1 != value_at_p2 && ! ( rdomain.null(value_at_p1) && rdomain.null(value_at_p2) ) ); const Surface_patch_index sp_index = rdomain.make_surface_index(value_at_p1, value_at_p2); const Index index = rdomain.index_from_surface_patch_index(sp_index); return Intersection(mid, index, 2); }

    // Else we must go on
    // Here we consider that p1(a) is the source point. Thus, we keep p1 and
    // change p2 if f(p1)!=f(p2).
    // That allows us to find the first intersection from a of [a,b] with
    // a surface.
    if ( value_at_p1 != value_at_mid &&
         ! ( r_domain_.null(value_at_p1) && r_domain_.null(value_at_mid) ) )
    {
      p2 = mid;
      value_at_p2 = value_at_mid;
    }
    else
    {
      p1 = mid;
      value_at_p1 = value_at_mid;
    }

    mid = midpoint(p1, p2);
    value_at_mid = r_domain_.function_(mid);
  }

The issue might be that the points are very large (recall, the SDF is limited to [-1,1] so it is unclear why vertices outside this box are created)

Point a(-93511495021678.5, 55518573485140.023, 75908032948603.141) Point b(-8429316601328.5234, 1704982732075.3398, 2885354838982.9434)

and logs (#define CGAL_MESH_3_VERBOSE, #define CGAL_MESH_3_VERY_VERBOSE) show that: insert_bad_facet(Facet { -2648552672213.6421 789240988277.20288 2762735743946.9658 0 2 0 1 -3966738333578.5615 -2058345041383.3494 3325365363490.9746 0 2 0 1 1075153782785.879 3093275916248.7319 5861172391787.7461 0 2 0 1 4th vertex in cell: -738115492532.64404 -4265990933791.1714 7024928602550.582 0 2 0 1 }

when interrupting the infinite loop, values of p1 and p2 are: p1(-67592449666003.141, 39125047539090.344, 53662738563152.422) p2(-67592449666003.148, 39125047539090.352, 53662738563152.43)

r_domain.squared_errorbound = 3.0000000000000001e-06

so the break condition (squared_distance(p1, p2) < rdomain.squared_errorbound) can never be fulfilled due to limited double precision

stefango74 commented 3 months ago

I found the original problem in my code - SDF distance was not computed correctly if outside and projected onto the [-1,1] bounding box, leading to these large values. Once fixed, the algorithm terminates nicely.

CGAL / cgal