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 #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.