Highly Degenerate Polyhedra from halfspace_intersection_3()

[SamSavitz]

First off, thank you for the wonderful package. I really like the header-only design, templating, and robustness. I have bolded potential issues below.

Issue Details

I am using what I think to be fairly reasonable input data to halfspace_intersection_3(). I decided to use exact rationals in my construction as double precision input fails with error:

terminate called after throwing an instance of 'CGAL::Assertion_exception'
  what():  CGAL ERROR: assertion violation!
Expr: boost::get<Point_3>(& *result) != nullptr
File: /usr/include/CGAL/Convex_hull_3/dual/halfspace_intersection_3.h
Line: 102
Explanation: halfspace_intersection_3: intersection is not a point

This is because the dual convex hull has a triangular facet with vertices that represent the planes:

 5.8637033051562737*x - 6*y + 2.9318516525781373*z - 1 == 0
-4.8989794855663549*x - 6*y - 2.4494897427831765*z - 1 == 0
 6.0000000000000027*x - 6*y + 3.0000000000000013*z - 1 == 0

According to CGAL::intersection, these planes do not intersect, despite {0, -1./6, 0} being an apparently obvious solution. When I use L = 48 (Warning: The runtime and memory use scale at least as L⁴.), it runs without error, but the output contains and uses obviously erroneous vertices such as

{-3.81607e+13, 7.63214e+13, -0.161617}

I'm now realizing that halfspace_intersection_with_constructions_3() seems to avoid this issue, but has the same degenerate output issues as I describe below .

Now, with exact constructions, I was getting back polyhedra that were so broken that I had to decompose them into polygon soup, repair, and then reconstruct them. (I tried stitching and duplicating vertices in boudary cycles, but, as my mesh is closed, those operations did not help.) I couldn't find an easier way to decompose a Polyhedron_3 into polygon soup, so I based my code on the .OFF format generator in Surface_mesh::write_off(). Perhaps there should be an easier or better documented way to convert a Polyhedron_3 to polygon soup. (Maybe this would be easier with a Surface_mesh, as that data structure is already index based.)

However, this was still insufficient, as simplify_polygon() in CGAL/Polygon_mesh_processing/repair_polygon_soup.h appears to fail when there are, in my case, three consecutive repeated vertices in the polygon. In my sample code fail.cpp, this first occurs for polygon_index == 61 in simplify_polygons_in_polygon_soup() where the input polygon == {65, 133, 137, 65, 65}. Replacing the body of simplify_polygon() with the code below seems to fix and shorten the function, although I do not know if it is a suitable drop-in replacement for all use cases.

Source Code

fail.cpp

simplify_polygon() Replacement

  const std::size_t ini_polygon_size = polygon.size();

  for (size_t i = 0; i < polygon.size(); i++) {
    size_t ni = (i + 1) % polygon.size();

    if (polygon[i] == polygon[ni] || traits.equal_3_object() (points[polygon[i]], points[polygon[ni]]))
      polygon.erase(polygon.begin() + i--);
  }

  const std::size_t removed_points_n = ini_polygon_size - polygon.size();
  return (removed_points_n != 0);

Stack trace with -O3 replaced with -ggdb

#0  __GI_raise (sig=sig@entry=6) at ../sysdeps/unix/sysv/linux/raise.c:51
#1  0x00007ffff718142a in __GI_abort () at abort.c:89
#2  0x00007ffff78810ad in __gnu_cxx::__verbose_terminate_handler() () from /usr/lib/x86_64-linux-gnu/libstdc++.so.6
#3  0x00007ffff787f066 in ?? () from /usr/lib/x86_64-linux-gnu/libstdc++.so.6
#4  0x00007ffff787f0b1 in std::terminate() () from /usr/lib/x86_64-linux-gnu/libstdc++.so.6
#5  0x00007ffff787f2c9 in __cxa_throw () from /usr/lib/x86_64-linux-gnu/libstdc++.so.6
#6  0x00005555555ea685 in boost::throw_exception<boost::bad_any_cast> (e=...) at /usr/include/boost/throw_exception.hpp:69
#7  0x00005555555c3ae5 in boost::any_cast<unsigned int> (operand=...) at /usr/include/boost/any.hpp:265
#8  0x00005555555c36df in CGAL::SNC_FM_decorator<CGAL::SNC_structure<CGAL::Epeck, CGAL::SNC_indexed_items, bool> >::determine_facet (this=0x7fffffffd710, e=..., MinimalEdge=std::vector of length 1, capacity 1 = {...}, FacetCycle=..., Edge_of=std::vector of length 5, capacity 5 = {...}) at /usr/include/CGAL/Nef_3/SNC_FM_decorator.h:416
#9  0x00005555555a0cec in CGAL::SNC_FM_decorator<CGAL::SNC_structure<CGAL::Epeck, CGAL::SNC_indexed_items, bool> >::create_facet_objects (this=0x7fffffffd710, plane_supporting_facet=..., start={obj = {px = 0x4, pn = {pi_ = 0x555556d9db60}}}, end={obj = {px = 0x4, pn = {pi_ = 0x555556d9db60}}}) at /usr/include/CGAL/Nef_3/SNC_FM_decorator.h:646
#10 0x0000555555587685 in CGAL::SNC_external_structure<CGAL::SNC_indexed_items, CGAL::SNC_structure<CGAL::Epeck, CGAL::SNC_indexed_items, bool> >::categorize_facet_cycles_and_create_facets (this=0x7fffffffd980) at /usr/include/CGAL/Nef_3/SNC_external_structure.h:1272
#11 0x00005555555782ee in CGAL::SNC_external_structure<CGAL::SNC_indexed_items, CGAL::SNC_structure<CGAL::Epeck, CGAL::SNC_indexed_items, bool> >::build_external_structure (this=0x7fffffffd980) at /usr/include/CGAL/Nef_3/SNC_external_structure.h:1335
#12 0x000055555556ee9d in CGAL::Nef_polyhedron_3<CGAL::Epeck, CGAL::SNC_indexed_items, bool>::build_external_structure (this=0x7fffffffdc00) at /usr/include/CGAL/Nef_polyhedron_3.h:349
#13 0x000055555556bcde in CGAL::Nef_polyhedron_3<CGAL::Epeck, CGAL::SNC_indexed_items, bool>::Nef_polyhedron_3<CGAL::Epeck, CGAL::Polyhedron_items_3, CGAL::HalfedgeDS_default, std::allocator<int> > (this=0x7fffffffdc00, P=...) at /usr/include/CGAL/Nef_polyhedron_3.h:604
#14 0x000055555556304f in main () at fail.cpp:45

Environment

Debian Stretch g++ -frounding-math -O3 -march=native -std=c++17 fail.cpp -o Fail -lgmpxx -lgmp -lmpfr CGAL version 5.0.2-3 installed from Debian's "sid" repository. Boost version 1.62.0.1 GMP version 2:6.2.0+dfsg-6 MPFR version 3.1.5-1

[MaelRL]

I have fixed the polygon soup repairing issue in #4884, thanks for the patch.

[lrineau]

I reopen, because https://github.com/CGAL/cgal/pull/4884 says it only "partially fixes" this issue.

[SamSavitz]

Thanks for the effort, @MaelRL and @lrineau. This has been my first contribution to an open source project, and you made it easy and pleasant.

As this was reopened, I figured I'd add my code for decomposing a Polyhedron_3 p into the polygon soup vs and fs, based on Surface_mesh::write_off():

vector <decltype(p)::Point> vs;
vs.reserve(num_vertices(p));

using VD = boost::graph_traits <decltype(p)> ::vertex_descriptor;
boost::container::flat_map <VD, size_t> reindex;
reindex.reserve(num_vertices(p));

const auto& vpm = get_const_property_map(CGAL::vertex_point, p);

size_t n = 0;
for (const auto& v:  vertices(p)) {
  vs.push_back(get(vpm, v));
  reindex[v] = n++;
}

vector <vector <size_t> > fs;
fs.reserve(num_faces(p));

transform(
  faces(p).begin(  ),
  faces(p).  end(  ),
  back_inserter (fs),
  [&, p, reindex] (const auto& f) {
    vector <size_t> face;
    face.reserve(f->facet_degree());

    const auto& vaf = vertices_around_face(halfedge(f, p), p);

    transform(
      vaf.begin    (    ),
      vaf.end      (    ),
      back_inserter(face),
      [&reindex] (const VD& vd) {
        return reindex.at(vd);
      });

    return face;
  });

Once again, I can't guarantee that this is suitable or even correct for all use cases.

[SamSavitz]

Furthermore, with double precision, halfspace_intersection_3()'s error message incorrectly reports that the planes mentioned above have no intersection. In fact, CGAL::intersection() returns a Line_3 in the form of the boost::variant type result_inter. This case is not being properly handled in build_primal_polyhedron(), and I do not know what the correct behavior should be.

[sloriot]

The graph combinatoric is correct but the embedding of the point in double are invalid due to rounding errors. There is not much we can do here. 3D snap rounding would be needed if you cannot use a Kernel with exact constructions. Repairing the output is also a possible solution.

[SamSavitz]

Well, one easy possibility is that the error message Explanation: halfspace_intersection_3: intersection is not a point could be split into cases where the intersection is null, a line, or a plane. My first interpretation of it was that the intersection was null, although I now see that it does not imply that.

[SamSavitz]

A further update: After increasing L to 96 and R to 2⁴⁰, Mathematica was still having trouble loading my .off file. I checked for missed duplicate vertices and didn't find any to double precison. But it turns out face #2122 is {6, 14429, 8105, 6, 2281, 12787, 3723}. This is the only face with remaining duplicated vertices. This is with exact constructions. I assume the ideal behavior of repair_polygon_soup() would be to split this polygon into a coplanar triangle and a square joined at 6 vertex. Is this correct? If so I might try patching it further.

Hopefully the mesh isn't broken in other ways. Mathematica had a lot more than just two faces missing.

It turns out that taking the vertices left over after removing duplicates and running convex_hull_3() on them produces a valid polyhedron that looks better than any of my previous results. It's a triangulation, unlike before, for better or worse. Perhaps this should be suggested somewhere in the documentation, or the halfspace intersection can be automatically checked for validity and run through convex_hull_3() before being reconstructed, if necessary. Is there any more efficient convex hull algorithm to use if we know that all the vertices are extreme points? It's somewhat surprising considering that I assume a convex hull algorithm is being used internally. Both halfspace_intersections_3() and halfspace_intersection_with_constructions_3() fail in a similar way, although the latter is actually slower because it uses more RAM and pushes into my swap memory. Also, I had to cast CGAL::ORIGIN to a Point_3 <K> only when using constructions, for some reason.

[sloriot]

Do you mean you have a case that is using exact constructions that produces a face having 2 vertices with an identical point associated (without rounding) ?

[SamSavitz]

I was being dumb and interpreting the first 6 as a vertex index when it was really saying that that face was a hexagon. I think the bug here lies with Mathematica. OpenCTM manages to triangulate it just fine. After triangulation, whether by convex_hull_3() or OpenCTM, Mathematica seems to handle it better, but it still reports the wrong number of vertices. I'm guessing that it doesn't like how close some of them are?

I guess the only remaining part of this issue is that halfspace_intersection_3() returns a polyhedron with duplicated vertices that requires decomposition to polygon soup and polyhedron reconstruction to fix. I feel like halfspace_intersection_3() should perhaps check for duplicated vertices and/or there should be a CGAL function which decomposes polyhedra into polygon soup. I posted a potential body for the latter above, but don't know how to properly write the function declaration with the right name, templates, traits, documentation etc.

[sloriot]

Side node: you can use triangulate_faces() instead of convex_hull_3() IIUC what you are doing.