Failure of intersection computation (five tetrahedra) #6568

stla · 2022-05-07T07:48:45Z

stla
May 7, 2022

Hello,

Here is the compound of five tetrahedra:

These are five tetrahedra centered at the origin. I didn't manage to compute their intersection. I can get the intersection of four of them, but not of the five ones. I use CGAL via R and the computation crashes the R session, without informative error messages.

I would be grateful if someone tries. Attached are the off files and below is the (untested) code.

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
#include <CGAL/Polygon_mesh_processing/corefinement.h>
#include <CGAL/Surface_mesh.h>
#include <fstream>
#include <iostream>

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Surface_mesh<K::Point_3> Mesh;

namespace PMP = CGAL::Polygon_mesh_processing;

int main(int argc, char* argv[]) {
  const std::string filename1 =
      (argc > 1) ? argv[1] : CGAL::data_file_path("th1.off");
  const std::string filename2 =
      (argc > 2) ? argv[2] : CGAL::data_file_path("th2.off");
  const std::string filename3 =
      (argc > 3) ? argv[3] : CGAL::data_file_path("th3.off");
  const std::string filename4 =
      (argc > 4) ? argv[4] : CGAL::data_file_path("th4.off");
  const std::string filename5 =
      (argc > 5) ? argv[5] : CGAL::data_file_path("th5.off");

  Mesh mesh1, mesh2, mesh3, mesh4, mesh5;
  if(!PMP::IO::read_polygon_mesh(filename1, mesh1)) {
    std::cerr << "Invalid input th1." << std::endl;
    return 1;
  }
  if(!PMP::IO::read_polygon_mesh(filename2, mesh2)) {
    std::cerr << "Invalid input th2." << std::endl;
    return 1;
  }
  if(!PMP::IO::read_polygon_mesh(filename3, mesh3)) {
    std::cerr << "Invalid input th3." << std::endl;
    return 1;
  }
  if(!PMP::IO::read_polygon_mesh(filename4, mesh4)) {
    std::cerr << "Invalid input th4." << std::endl;
    return 1;
  }
  if(!PMP::IO::read_polygon_mesh(filename5, mesh5)) {
    std::cerr << "Invalid input th5." << std::endl;
    return 1;
  }

  Mesh inter12;
  if(PMP::corefine_and_compute_intersection(mesh1, mesh2, inter12)) {
    std::cout << "Intersection th1-th2 successfully computed.\n";
    Mesh inter34;
    if(PMP::corefine_and_compute_intersection(mesh3, mesh4, inter34)) {
      std::cout << "Intersection th3-th4 successfully computed.\n";
      Mesh inter1234;
      if(PMP::corefine_and_compute_intersection(inter12, inter34, inter1234)) {
        std::cout << "Intersection th1-th2-th3-th4 successfully computed.\n";
        Mesh inter12345;
        if(PMP::corefine_and_compute_intersection(inter1234, mesh5,
                                                  inter12345)) {
          std::cout << "Final intersection successfully computed.\n";
          CGAL::IO::write_polygon_mesh("inter_tetrahedra.off", inter12345,
                                       CGAL::parameters::stream_precision(17));
          return 0;
        } else {
          std::cout << "Final intersection failed.\n";
          return 1;
        }
      } else {
        std::cout << "Intersection th1-th2-th3-th4 failed.\n";
        return 1;
      }
    } else {
      std::cout << "Intersection th3-th4 failed.\n";
      return 1;
    }
  } else {
    std::cout << "Intersection th1-th2 failed.\n";
    return 1;
  }
}

If you try and this works, that would perhaps mean there's a problem in my R/C++ code. I also tried with nef polyhedras but this was worst.

ths_offFiles.zip

stla · 2022-05-07T07:55:34Z

stla
May 7, 2022
Author

Wait... writing the off files has rounded the vertices up to 4 digits, and this works with the rounded vertices! I get a pretty polyhedra for the intersection. Do you know why the rounding solves the issue?

0 replies

stla · 2022-05-07T07:57:17Z

stla
May 7, 2022
Author

Also, that's nice if the rounding is a solution, but I'd like to prevent the crash.

0 replies

stla · 2022-05-07T08:06:33Z

stla
May 7, 2022
Author

Attached are the off files with the unrounded vertices.
ths_unrounded.zip

0 replies

pentacular · 2022-05-07T10:09:27Z

pentacular
May 7, 2022

You're using an Exact_predicates_inexact_constructions_kernel

corefine_and_compute_intersection may produce self-intersecting results with this type of kernel.

Try using an Exact_predicates_exact_constructions_kernel if you want to reprocess the results as above.

2 replies

stla May 7, 2022
Author

Thank you, I will try soon. However, I would like to leave the choice to the R user: exact or inexact. So I'm still wondering whether there is a way to prevent the crash.

stla May 7, 2022
Author

I did implement with the exact kernel. Same result: crash.

afabri · 2022-05-08T13:21:26Z

afabri
May 8, 2022
Collaborator

Can you please give more context. I tried it with the master branch as well as with 5.4 with Visual C++ and boost 1_66_0 and I cannot reproduce your problem yet.

1 reply

stla May 9, 2022
Author

Thank you. It's a bit complicated to give more context since a priori you are not a R user, but let's go. It has always been possible to include some C code in R. Now a guy has made a package, Rcpp, which is an amazing piece of work allowing to include C++ in R. The same guy also maintains the R packages RcppEigen, which provides the latest version of Eigen, and BH, which provides the latest version of Boost. And someone else recently did the RcppCGAL package (currently providing CGAL 5.4). So I'm using all these packages. Maybe there's something bad in my code, I'm not very knowledgeable in C++, but it's strange that my package does not crash the R session when the vertex coordinates are rounded, or when I intersect only two tetrahedra. With Rcpp, the C++ code is compiled to a DLL in Windows with g++.

stla · 2022-05-09T09:03:14Z

stla
May 9, 2022
Author

I tried with a newer version of g++ and now this crashes even for two tetrahedra with rounded vertices. I'm trying on Linux now, via a Github action. That does not smell good: the action is running for 30 minutes, while it took 25 minutes before.

0 replies

stla · 2022-05-09T10:47:18Z

stla
May 9, 2022
Author

I tested with Windows on another machine, including a DEBUG option, and I get this error:

> CGAL ERROR: assertion violation!

 Expr: p1p2p3_2d!=COLLINEAR || p2p1p4_2d!=COLLINEAR

 File: C:/Users/USERqCzoAgSijr/AppData/Local/R/win-library/4.2/RcppCGAL/include/CGAL/Polygon_mesh_processing/orientation.h
> Line: 117

0 replies

afabri · 2022-05-09T10:57:06Z

afabri
May 9, 2022
Collaborator

I am sorry but you must be more precise. You finally do it with which kernel, and with which data set. What does "with Windows" mean, that is in case it is VC++ which version? What do you mean with "including a DEBUG option"? Also what does "4.2" in the line starting with "File: C:/" refer to ? Is that the release number of CGAL you use?

1 reply

stla May 9, 2022
Author

Sorry. For my latest tests, I tried with the Exact_predicates_inexact_constructions kernel. I'm using the vertices of the five tetrahedra, with unrounded coordinates. The Rcpp code does not compile with VC++ nor cygwin. We just compile it with g++. The DEBUG option is the C++ flag -UNDEBUG. When using this flag, the R session does not crash, it prints the above error messages instead. "4.2" is the version of R. The CGAL version in RcppCGAL is 5.4.

stla · 2022-05-09T12:07:31Z

stla
May 9, 2022
Author

Wait... it has worked on another Windows machine using the vertex coordinates rounded up to 11 digits. And I see in my latest OFF files that the coordinates are rounded up to 11 digits. I'm going to prepare a code with totally unrounded coordinates this time, hoping someone will try it.

0 replies

stla · 2022-05-09T12:35:06Z

stla
May 9, 2022
Author

This time this is maximally unrounded. Please anyone, could you try?

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
#include <CGAL/Polygon_mesh_processing/corefinement.h>
#include <CGAL/Surface_mesh.h>
#include <math.h>
#include <fstream>
#include <iostream>

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Surface_mesh<K::Point_3> Mesh;

namespace PMP = CGAL::Polygon_mesh_processing;

int main(int argc, char* argv[]) {
  double phi = (1.0 + sqrt(5.0)) / 2.0;
  double a = 1.0 / sqrt(3.0);
  double b = a / phi;
  double c = a * phi;

  Mesh mesh1, mesh2, mesh3, mesh4, mesh5;

  Mesh::Vertex_index v0 = mesh1.add_vertex(Point_3(0, b, c));
  Mesh::Vertex_index v1 = mesh1.add_vertex(Point_3(b, -c, 0));
  Mesh::Vertex_index v2 = mesh1.add_vertex(Point_3(a, a, -a));
  Mesh::Vertex_index v3 = mesh1.add_vertex(Point_3(-c, 0, b));
  Mesh::Vertex_index v4 = mesh2.add_vertex(Point_3(a, -a, -a));
  Mesh::Vertex_index v5 = mesh2.add_vertex(Point_3(a, a, a));
  Mesh::Vertex_index v6 = mesh2.add_vertex(Point_3(-a, -a, a));
  Mesh::Vertex_index v7 = mesh2.add_vertex(Point_3(-a, a, -a));
  Mesh::Vertex_index v8 = mesh3.add_vertex(Point_3(c, 0, b));
  Mesh::Vertex_index v9 = mesh3.add_vertex(Point_3(-a, a, a));
  Mesh::Vertex_index v10 = mesh3.add_vertex(Point_3(-b, -c, 0));
  Mesh::Vertex_index v11 = mesh3.add_vertex(Point_3(0, b, -c));
  Mesh::Vertex_index v12 = mesh4.add_vertex(Point_3(a, -a, a));
  Mesh::Vertex_index v13 = mesh4.add_vertex(Point_3(b, c, 0));
  Mesh::Vertex_index v14 = mesh4.add_vertex(Point_3(-c, 0, b));
  Mesh::Vertex_index v15 = mesh4.add_vertex(Point_3(0, -b, -c));
  Mesh::Vertex_index v16 = mesh5.add_vertex(Point_3(-a, -a, -a));
  Mesh::Vertex_index v17 = mesh5.add_vertex(Point_3(-b, c, 0));
  Mesh::Vertex_index v18 = mesh5.add_vertex(Point_3(c, 0, -b));
  Mesh::Vertex_index v19 = mesh5.add_vertex(Point_3(0, -b, c));

  mesh1.add_face(v0, v1, v2);
  mesh1.add_face(v2, v1, v3);
  mesh1.add_face(v3, v1, v0);
  mesh1.add_face(v0, v2, v3);
  mesh2.add_face(v4, v5, v6);
  mesh2.add_face(v6, v5, v7);
  mesh2.add_face(v7, v5, v4);
  mesh2.add_face(v4, v6, v7);
  mesh3.add_face(v8, v9, v10);
  mesh3.add_face(v10, v9, v11);
  mesh3.add_face(v11, v9, v8);
  mesh3.add_face(v8, v10, v11);
  mesh4.add_face(v12, v13, v14);
  mesh4.add_face(v14, v13, v15);
  mesh4.add_face(v15, v13, v12);
  mesh4.add_face(v12, v14, v15);
  mesh5.add_face(v16, v17, v18);
  mesh5.add_face(v18, v17, v19);
  mesh5.add_face(v19, v17, v16);
  mesh5.add_face(v16, v18, v19);

  Mesh inter12;
  if(PMP::corefine_and_compute_intersection(mesh1, mesh2, inter12)) {
    std::cout << "Intersection th1-th2 successfully computed.\n";
    if(PMP::does_self_intersect(inter12)) {
      std::cout << "Intersection th1-th2 self-intersects.\n";
      return 0;
    }
    if(!PMP::does_bound_a_volume(inter12)) {
      std::cout << "Intersection th1-th2 does not bound a volume.\n";
      return 0;
    }
    Mesh inter34;
    if(PMP::corefine_and_compute_intersection(mesh3, mesh4, inter34)) {
      std::cout << "Intersection th3-th4 successfully computed.\n";
      if(PMP::does_self_intersect(inter34)) {
        std::cout << "Intersection th3-th4 self-intersects.\n";
        return 0;
      }
      if(!PMP::does_bound_a_volume(inter34)) {
        std::cout << "Intersection th3-th4 does not bound a volume.\n";
        return 0;
      }
      Mesh inter1234;
      if(PMP::corefine_and_compute_intersection(inter12, inter34, inter1234)) {
        std::cout << "Intersection th1-th2-th3-th4 successfully computed.\n";
        if(PMP::does_self_intersect(inter1234)) {
          std::cout << "Intersection th1-th2-th3-th4 self-intersects.\n";
          return 0;
        }
        if(!PMP::does_bound_a_volume(inter1234)) {
          std::cout
              << "Intersection th1-th2-th3-th4 does not bound a volume.\n";
          return 0;
        }
        Mesh inter12345;
        if(PMP::corefine_and_compute_intersection(inter1234, mesh5,
                                                  inter12345)) {
          std::cout << "Final intersection successfully computed.\n";
          CGAL::IO::write_polygon_mesh("inter_tetrahedra.off", inter12345,
                                       CGAL::parameters::stream_precision(17));
          return 0;
        } else {
          std::cout << "Final intersection failed.\n";
          return 1;
        }
      } else {
        std::cout << "Intersection th1-th2-th3-th4 failed.\n";
        return 1;
      }
    } else {
      std::cout << "Intersection th3-th4 failed.\n";
      return 1;
    }
  } else {
    std::cout << "Intersection th1-th2 failed.\n";
    return 1;
  }
}

17 replies

sloriot May 10, 2022
Maintainer

point me to your wrapper code in github and I'll try to have a look.

stla May 10, 2022
Author

Thanks. Maybe the exception is caught because I set the compilation option -UNDEBUG (this option prevents the crashes)? I unset it and I ran gdb one hour ago and it is still running.
My code is here, function Intersection2.

sloriot May 13, 2022
Maintainer

Line 285, you are writing the output in meshes[i] which will be first mesh at the next loop. You are not checking that it is free from self-intersection (which is not needed with a Kernel with exact constructions). For iterative Boolean operations, I would recommend to use that example. Note that you may observe a slow down after a few iteration due to the increase in the size of the multiprecision numbers. An alternative is to clean the mesh when it is self-intersecting with for example CGAL::Polygon_mesh_processing::experimental::remove_self_intersections().

stla May 13, 2022
Author

Thanks. Ok, so you mean it is possible that the computed intersection self-intersects while the two input meshes do not. I didn't know that, thanks.

sloriot May 13, 2022
Maintainer

when using an inexact kernel, you have some rounding going on and self-intersections can happen.

afabri · 2022-05-10T10:33:51Z

afabri
May 10, 2022
Collaborator

And please put the code on gist.github.com

0 replies

stla · 2022-05-12T14:18:52Z

stla
May 12, 2022
Author

Failure confirmed: darikg/seagullmesh#1

2 replies

sloriot May 13, 2022
Maintainer

Converted the code into C++ and it does not crash:
https://gist.github.com/sloriot/95c33925653f120b086a987cdb295371

However, due to the rounding the output is self-intersecting.

stla May 13, 2022
Author

It works now !! I'm lost, I don't remember I changed something...

stla · 2022-05-14T13:12:58Z

stla
May 14, 2022
Author

Blog post: https://laustep.github.io/stlahblog/posts/BooleanOpsOnMeshes.html

0 replies

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Failure of intersection computation (five tetrahedra) #6568

{{title}}

Replies: 13 comments 23 replies

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

Select a reply

Failure of intersection computation (five tetrahedra) #6568

stla May 7, 2022

Replies: 13 comments · 23 replies

stla May 7, 2022 Author

stla May 7, 2022 Author

stla May 7, 2022 Author

pentacular May 7, 2022

stla May 7, 2022 Author

stla May 7, 2022 Author

afabri May 8, 2022 Collaborator

stla May 9, 2022 Author

stla May 9, 2022 Author

stla May 9, 2022 Author

afabri May 9, 2022 Collaborator

stla May 9, 2022 Author

stla May 9, 2022 Author

stla May 9, 2022 Author

sloriot May 10, 2022 Maintainer

stla May 10, 2022 Author

sloriot May 13, 2022 Maintainer

stla May 13, 2022 Author

sloriot May 13, 2022 Maintainer

afabri May 10, 2022 Collaborator

stla May 12, 2022 Author

sloriot May 13, 2022 Maintainer

stla May 13, 2022 Author

stla May 14, 2022 Author

stla
May 7, 2022

Replies: 13 comments 23 replies

stla
May 7, 2022
Author

stla
May 7, 2022
Author

stla
May 7, 2022
Author

pentacular
May 7, 2022

stla May 7, 2022
Author

stla May 7, 2022
Author

afabri
May 8, 2022
Collaborator

stla May 9, 2022
Author

stla
May 9, 2022
Author

stla
May 9, 2022
Author

afabri
May 9, 2022
Collaborator

stla May 9, 2022
Author

stla
May 9, 2022
Author

stla
May 9, 2022
Author

sloriot May 10, 2022
Maintainer

stla May 10, 2022
Author

sloriot May 13, 2022
Maintainer

stla May 13, 2022
Author

sloriot May 13, 2022
Maintainer

afabri
May 10, 2022
Collaborator

stla
May 12, 2022
Author

sloriot May 13, 2022
Maintainer

stla May 13, 2022
Author

stla
May 14, 2022
Author