CGAL / cgal

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

CGAL::Euler::add_face() creates a non-manifold edge #4919

Open sgiraudot opened 4 years ago

sgiraudot commented 4 years ago

Issue Details

CGAL::Euler_add_face() is supposed to return a null_face if the face could not be created (typically when a non-manifold edge would be created).

For some reason, this particular case makes CGAL::Euler::add_face() create a non-manifold edge without protesting:

image

Adding face 6 is obviously wrong (edge 0-4 has 3 neighbors) but it does not trigger any assertion. As a result, trying to load this mesh in a CGAL::Surface_mesh works without triggering any error. Calling mesh::is_valid() returns true.

(Of course in that case, there are also non-manifold vertices, but these can be stored in a CGAL::Surface_mesh and are supposed to be handled separately, so I'm less surprised).

When calling CGAL::PMP::is_polygon_soup_a_polygon_mesh() on the created CGAL::Surface_mesh, it does detect the non-manifold edge and thus does return false, which is at least correct.

I'm mainly posting this issue for keeping track of it, but I'm still investigating.

Source Code

I've reduce the bug to the mesh corresponding to the screenshot above:

OFF
10 7 0
0 0 0
0 1 -0.5
0 1 -1
0 -1 0
0 0 1
0 1 0.5
0 1 0
0.75 0 0.5
0.5 -0.5 0
0 0.25 0.5
3 0 1 2
3 0 3 4
3 5 6 0
3 0 7 8
3 4 7 0
3 3 0 8
3 9 4 0

The following source code either reads this file and loads a CGAL::Surface_mesh, or it constructs the same mesh (without geometry) using BGL + Euler operations:

#include <iostream>
#include <fstream>

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
#include <CGAL/IO/OFF_reader.h>

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

struct Face_to_index_vector
{
  typedef Mesh::Face_index argument_type;
  typedef std::vector<std::size_t> return_type;

  const Mesh& mesh;

  Face_to_index_vector (const Mesh& mesh) : mesh(mesh) { }

  return_type operator() (const argument_type& fi) const
  {
    return_type out;
    for(typename Mesh::Vertex_index v : CGAL::vertices_around_face(mesh.halfedge(fi),mesh))
      out.push_back (std::size_t(v));
    return out;
  }
};

int main(int argc, char** argv)
{
  if (argc > 1)
  {
    std::ifstream file(argv[1]);

    Mesh mesh;
    if (!CGAL::read_off(file, mesh)) // This should fail but doesn't
      std::cerr << "Can't read" << std::endl;

    if (file.fail()) // This should then be true (but it's not)
      std::cerr << "File reading failure" << std::endl;
    if (mesh.is_empty()) // Then it should be empty (but it's not)
      std::cerr << "Empty" << std::endl;
    if (!mesh.is_valid()) // Or at least it should be invalid (but it's not)
      std::cerr << "Invalid" << std::endl;

    // This DOES fail, the mesh is not a polygon mesh (!)
    Face_to_index_vector f2i(mesh);
    if (!CGAL::Polygon_mesh_processing::is_polygon_soup_a_polygon_mesh
        (boost::make_iterator_range (boost::make_transform_iterator(mesh.faces().begin(), f2i),
                                     boost::make_transform_iterator(mesh.faces().end(), f2i))))
      std::cerr << "Not a mesh!" << std::endl;
  }
  else
  {
    Mesh mesh;
    auto v0 = CGAL::add_vertex(mesh);
    auto v1 = CGAL::add_vertex(mesh);
    auto v2 = CGAL::add_vertex(mesh);
    auto v3 = CGAL::add_vertex(mesh);
    auto v4 = CGAL::add_vertex(mesh);
    auto v5 = CGAL::add_vertex(mesh);
    auto v6 = CGAL::add_vertex(mesh);
    auto v7 = CGAL::add_vertex(mesh);
    auto v8 = CGAL::add_vertex(mesh);
    auto v9 = CGAL::add_vertex(mesh);

    auto f0 = CGAL::Euler::add_face (CGAL::make_array(v0, v1, v2), mesh);
    CGAL_assertion (f0 != boost::graph_traits<Mesh>::null_face());
    auto f1 = CGAL::Euler::add_face (CGAL::make_array(v0, v3, v4), mesh);
    CGAL_assertion (f1 != boost::graph_traits<Mesh>::null_face());
    auto f2 = CGAL::Euler::add_face (CGAL::make_array(v5, v6, v0), mesh);
    CGAL_assertion (f2 != boost::graph_traits<Mesh>::null_face());
    auto f3 = CGAL::Euler::add_face (CGAL::make_array(v0, v7, v8), mesh);
    CGAL_assertion (f3 != boost::graph_traits<Mesh>::null_face());
    auto f4 = CGAL::Euler::add_face (CGAL::make_array(v4, v7, v0), mesh);
    CGAL_assertion (f4 != boost::graph_traits<Mesh>::null_face());
    auto f5 = CGAL::Euler::add_face (CGAL::make_array(v3, v0, v8), mesh);
    CGAL_assertion (f5 != boost::graph_traits<Mesh>::null_face());

    // This should obviously fail as it creates a non-manifold edge 0-4 (with 3 neighbor faces f1, f4 and f6)
    auto f6 = CGAL::Euler::add_face (CGAL::make_array(v9, v4, v0), mesh);
    CGAL_assertion (f6 == boost::graph_traits<Mesh>::null_face());
  }

  return EXIT_SUCCESS;
}

Environment

sgiraudot commented 4 years ago

The bug is quite identified now: when closing the "umbrella" adding f5, the HDS is broken, v4-v0 is no longer considered an halfedge. Closing this umbrella actually makes the non-manifoldness of v0 "unsolvable" by adding new faces and it should probably return a null face at this point.

Code with updated assertions:

#include <iostream>
#include <fstream>

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
#include <CGAL/IO/OFF.h>

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

struct Face_to_index_vector
{
  typedef Mesh::Face_index argument_type;
  typedef std::vector<std::size_t> return_type;

  const Mesh& mesh;

  Face_to_index_vector (const Mesh& mesh) : mesh(mesh) { }

  return_type operator() (const argument_type& fi) const
  {
    return_type out;
    for(typename Mesh::Vertex_index v : CGAL::vertices_around_face(mesh.halfedge(fi),mesh))
      out.push_back (std::size_t(v));
    return out;
  }
};

int main(int argc, char** argv)
{
  if (argc > 1)
  {
    std::ifstream file(argv[1]);

    Mesh mesh;
    if (!CGAL::IO::read_OFF(file, mesh))
      std::cerr << "Can't read" << std::endl;
    if (file.fail()) // This should fail but doesn't
      std::cerr << "File reading failure" << std::endl;
    if (mesh.is_empty()) // Then it should be empty (but it's not)
      std::cerr << "Empty" << std::endl;
    if (!mesh.is_valid()) // Or at least it should be invalid (but it's not)
      std::cerr << "Invalid" << std::endl;

    // This DOES fail, the mesh is not a polygon mesh (!)
    Face_to_index_vector f2i(mesh);
    if (!CGAL::Polygon_mesh_processing::is_polygon_soup_a_polygon_mesh
        (boost::make_iterator_range (boost::make_transform_iterator(mesh.faces().begin(), f2i),
                                     boost::make_transform_iterator(mesh.faces().end(), f2i))))
      std::cerr << "Not a mesh!" << std::endl;
  }
  else
  {
    Mesh mesh;
    auto v0 = CGAL::add_vertex(mesh);
    auto v1 = CGAL::add_vertex(mesh);
    auto v2 = CGAL::add_vertex(mesh);
    auto v3 = CGAL::add_vertex(mesh);
    auto v4 = CGAL::add_vertex(mesh);
    auto v5 = CGAL::add_vertex(mesh);
    auto v6 = CGAL::add_vertex(mesh);
    auto v7 = CGAL::add_vertex(mesh);
    auto v8 = CGAL::add_vertex(mesh);
    auto v9 = CGAL::add_vertex(mesh);

    auto f0 = CGAL::Euler::add_face (CGAL::make_array(v0, v1, v2), mesh);
    CGAL_assertion (f0 != boost::graph_traits<Mesh>::null_face());
    auto f1 = CGAL::Euler::add_face (CGAL::make_array(v0, v3, v4), mesh);
    CGAL_assertion (f1 != boost::graph_traits<Mesh>::null_face());

    CGAL_assertion (halfedge(v4, v0, mesh).second);
    CGAL_assertion (halfedge(v0, v4, mesh).second);
    CGAL_assertion (!CGAL::is_border(halfedge(v4, v0, mesh).first, mesh));
    CGAL_assertion (CGAL::is_border(halfedge(v0, v4, mesh).first, mesh));

    auto f2 = CGAL::Euler::add_face (CGAL::make_array(v5, v6, v0), mesh);
    CGAL_assertion (f2 != boost::graph_traits<Mesh>::null_face());
    auto f3 = CGAL::Euler::add_face (CGAL::make_array(v0, v7, v8), mesh);
    CGAL_assertion (f3 != boost::graph_traits<Mesh>::null_face());
    auto f4 = CGAL::Euler::add_face (CGAL::make_array(v4, v7, v0), mesh);
    CGAL_assertion (f4 != boost::graph_traits<Mesh>::null_face());

    CGAL_assertion (halfedge(v4, v0, mesh).second);
    CGAL_assertion (halfedge(v0, v4, mesh).second);
    CGAL_assertion (!CGAL::is_border(halfedge(v4, v0, mesh).first, mesh));
    CGAL_assertion (!CGAL::is_border(halfedge(v0, v4, mesh).first, mesh));

    auto f5 = CGAL::Euler::add_face (CGAL::make_array(v3, v0, v8), mesh);
    CGAL_assertion (f5 != boost::graph_traits<Mesh>::null_face());

    // Already broken here -> v4 v0 is no longer a halfedge
    CGAL_assertion (halfedge(v4, v0, mesh).second);
    CGAL_assertion (halfedge(v0, v4, mesh).second);
    CGAL_assertion (!CGAL::is_border(halfedge(v4, v0, mesh).first, mesh));
    CGAL_assertion (!CGAL::is_border(halfedge(v0, v4, mesh).first, mesh));

    // This should obviously fail as it creates a non-manifold edge 0-4 (with 3 neighbor faces f1, f4 and f6)
    auto f6 = CGAL::Euler::add_face (CGAL::make_array(v9, v4, v0), mesh);
    CGAL_assertion (f6 == boost::graph_traits<Mesh>::null_face());
  }

  return EXIT_SUCCESS;
}
sloriot commented 4 years ago

The issue comes from a non-manifold vertex that prevent the existing edge to be found (before we are in the wrong umbrella). I don't see any other solution but to forbid non-manifold vertices. @MaelRL we need to discuss that when you're back.