search

CGAL / cgal

The public CGAL repository, see the README below
https://github.com/CGAL/cgal#readme
Other
4.9k stars 1.38k forks source link

Cannot clip mesh with cavity: Assertion_exception on s != LEFT_TURN Triangulation_2.h:931 #7493

Closed mildred closed 1 year ago

mildred commented 1 year ago

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

Issue Details

I have a complex geometry using Surface_mesh. It is constructed by boolean operations and it is basically a solid rectangle which is hollowed inside (in rectangle form), only a thin shell is around. Within the inner rectangle there is another rectangle added that makes a bump.

2023-06-03_21-50

mesh gzipped: dbg_convex_step_0.obj.gz

I'm trying to split this geometry to form multiple convex meshes. I'm using clip() with clip_volume(true). I'm clipping with a plane normal to the Z axis (0,0,1) at the offset z=-128. This should open up one face of the shell as the z=-128 coordinate is coplanar with the inner substracted space.

The exception message:

Clip #0 mesh 1 / 1
get first half-edge
get second half-edge
get opposite points (-512 512 -128), (-128 -16 -64)
get face normal 0 0 1
dot product: 64
clip current mesh
terminate called after throwing an instance of 'CGAL::Assertion_exception'
  what():  CGAL ERROR: assertion violation!
Expr: s != LEFT_TURN
File: /usr/include/CGAL/Triangulation_2.h
Line: 931

Program received signal SIGABRT, Aborted.

basic stack trace:

CGAL_triangulation_assertion( s != LEFT_TURN ); fails within the is_valid() function. In Visitor.h:1406 is_valid() is called on the variable cdt, and this is this variable that is the cause of the crash. When this variable is defined, there is a strange TODO comment:

///TODO use a positive normal and remove all workaround to guarantee that triangulation of coplanar patches are compatible
      CDT_traits traits(typename EK::Construct_normal_3()(p,q,r));
      CDT cdt(traits);

Source Code

To give you an idea of the code:

  bool clip_mesh(Mesh &current, std::list<Mesh> &queue, std::list<Mesh> &res) {
    using namespace CGAL::Polygon_mesh_processing;
    using namespace CGAL::parameters;
    using namespace CGAL;
    using Vector_3 = Kernel::Vector_3;
    using Plane_3 = Kernel::Plane_3;

    if (!is_triangle_mesh(current)) {
      cerr << "Mesh is not triangulated" << endl;
      return false;
    }
    if (!is_valid_polygon_mesh(current)) {
      cerr << "Mesh is not valid polygon" << endl;
      return false;
    }

    bool convex = true;

    // iterate over all faces
    for(edge_descriptor edge : current.edges()) {
      cerr << "get first half-edge" << endl;
      // get the first half edge and 3 vertices around
      auto edge1 = current.halfedge(edge);
      auto v1p = current.target(current.prev(edge1));
      auto v1t = current.target(edge1);
      auto v1n = current.target(current.next(edge1));

      cerr << "get second half-edge" << endl;
      // get the second half edge and 3 vertices around
      auto edge2 = current.opposite(edge1);
      auto v2p = current.target(current.prev(edge2));
      auto v2t = current.target(edge2);
      auto v2n = current.target(current.next(edge2));

      // get opposite points for each of the two faces
      // these points do not belongs to the edge
      auto p1 = current.point((v1p == v2t) ? v1n : v1p);
      auto p2 = current.point((v2p == v1t) ? v2n : v2p);
      cerr << "get opposite points (" << p1 << "), (" << p2 << ")" << endl;

      // compute the normal for the reference face
      auto norm1 = compute_face_normal(current.face(edge1), current);
      cerr << "get face normal " << norm1 << endl;

      // if the scalar product is positive, the angle is concave
      // <https://stackoverflow.com/a/40019587>
      auto product = scalar_product(Vector_3(p1, p2), norm1);
      cerr << "dot product: " << product << endl;
      if( product > 0 ) {
        Plane_3 clip_plane(p1, norm1);
        Plane_3 rclip_plane(p1, -norm1);

        convex = false;
        Mesh complement(current);
        cerr << "clip current mesh" << endl;
        clip(current, clip_plane, clip_volume(true));
        cerr << "clip complement mesh" << endl;
        clip(complement, rclip_plane, clip_volume(true));
        cerr << "push back 2 meshes" << endl;
        queue.push_back(complement);
        queue.push_back(current);
        // If the face iterator can work with a mesh changing, perhaps we can
        // continue to iterate over it but in the meantime put both meshes for
        // reprocess
        return true;
      }
    }

    if (convex) {
      // No concave surface got split, the mesh is convex
      res.push_back(current);
    }

    return true;
  }

source code: https://github.com/mildred/t3d2map/commit/53bc8c6815278bb098abc468f8ea7a9397114502

command-line: make && rm -f dbg_*.obj && ./t3d2map --debug-mesh ./examples/test.t3d

complete source code is one file, not that long. First part is constructing the geometry with CSG then

Environment

mildred commented 1 year ago

I'm trying to replace clip with a standard corefine_and_compute_difference with a large bounding box clipped by my clip plane and I got an assertion error too. Could it come from the geometry that is too unusual?

terminate called after throwing an instance of 'CGAL::Assertion_exception'
  what():  CGAL ERROR: assertion violation!
Expr: insert_ok || mesh_to_vertex_to_node_id[&tm][target(h,tm)]==node_id
File: /usr/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Visitor.h
Line: 768
Error: signal 6:
sloriot commented 1 year ago

Just looking at the mesh, I can see that it is not correctly oriented (nested volume should be counterclockwise oriented). If you call the function PMP::does_bound_a_volume() it should return false. Calling the function PMP::orient_to_bound_a_volume() should fix it. Please let us know if it is enough to fix the issue.

mildred commented 1 year ago

The normals from the cavity are oriented to the inside, which is away from the solid volume itself, which seems completely logical to me. That's how the boolean operations oriented them anyway. (By the way that's not a nested volume but a cavity).

After calling orient_to_bound_a_volume, the normals do not seems flipped (see dbg_convex_step_0.obj.gz), and I added an extra check with does_bound_a_volume before the clip but I got the same exception.

The volume contains a cavity, perhaps clip is not able to deal with that. But how different is a cavity from let's say a glass shaped convex volume when clipping?

sloriot commented 1 year ago

Sorry the input is indeed valid. However, you should not use using Kernel = CGAL::Simple_cartesian<double>; but using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel. If you know that the plane is axis aligned I would recommend the second approach you tried with a bbox. With the plane, I noticed that due to rounding issues, some coordinates of the internally created box have -128.00000000000003 coordinates which lead to an output which is not what you would expect.

mildred commented 1 year ago

Using CGAL::Exact_predicates_inexact_constructions_kernel instead of CGAL::Simple_cartesian<double> made it work. Thank you. I didn't know the kernel could have such dramatic consequences.

As for the rounding, this is to be expected but perhaps I can round all coordinates to some limited precision so I can get back neatly arranged vertices. And for the cutting plane, I have no guarantee that it should match an axis, so the bounding box creation involves some cutting in any case along a possibly non orthogonal plane.

mildred commented 1 year ago

What is the reason of the assertion failure, can I expect assertion failures due to geometry or am I safe with the new kernel?

edit: I think I got my anser here: https://www.cgal.org/exact.html

Thank you very much. I'll let you close the issue.