CGAL users discussion list

[Benjamin Kehlet <>]

Dear CGAL team

Thank you very much for your work with the Surface_mesh and Polygon_mesh processing_packages, and in particular the newly added implementation of isotropic remeshing. It is much appreciated!

When playing with the remesher the output is occasionally problematic. Attached is a simple program which reproduces the problem. It reads in a polyhedron from file, marks some edges as protected constraints, invokes the remesher and outputs to off format.

When running the program with the attached cylinder.off as input, the generated file (remeshed.off) contains:

Line 320: 3 126 124 122

Line 530: 3 122 124 126

Line 532: 3 281 122 126

Line 827: 3 294 126 122

So the faces in line 320 and 530 seems to be the same but with opposite orientation. The faces in line 320 and 532 share a halfedge from v122 to v126 and the faces in line 530 and 827 shared a halfedge from v126 to v122.

I'm on debian testing with a debug build of CGAL 4.8. I'm building with gcc 5.3.1 (with -std=gnu++11).

I hope someone can take a look an comment if I'm doing something wrong, the input can not be handled by the function or if this is a problem in CGAL.

Thanks in advance!

Best regards

Benjamin Kehlet

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/boost/graph/graph_traits_Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/remesh.h>
#include <CGAL/Polygon_mesh_processing/border.h>
#include <boost/function_output_iterator.hpp>
#include <fstream>
#include <vector>

typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::FT FT;
typedef Kernel::Point_3 Point;
typedef Kernel::Vector_3 Vector;
typedef CGAL::Surface_mesh<Point> Surface_mesh;
typedef Surface_mesh::Vertex_index Vertex_index;
typedef Surface_mesh::Face_index Face_index;
typedef Surface_mesh::Edge_index Edge_index;
typedef Surface_mesh::Halfedge_index Halfedge_index;
typedef boost::graph_traits<Surface_mesh>::face_descriptor face_descriptor;
typedef boost::graph_traits<Surface_mesh>::edge_descriptor edge_descriptor;


int main(int argc, char** argv)
{
  namespace Parameters = CGAL::Polygon_mesh_processing::parameters;

  Surface_mesh m;

  {
    std::ifstream infile(argv[1]);
    infile >> m;
   }

  std::cout << "Num vertices: " << m.number_of_vertices() << std::endl;
  std::cout << "Num faces:    " << m.number_of_faces() << std::endl;

  //----------------  Protect egdes sharper than 60 degrees ------------------
  const double cos_tolerance = .5;

  Surface_mesh::Property_map<edge_descriptor, bool> edge_constrained_map =
    m.add_property_map<edge_descriptor, bool>("e:is_constrained", false).first;

  std::size_t num_protected = 0;
  double longest_protected = 0.;
  Surface_mesh::Property_map<face_descriptor, Vector> fnormals =
    m.add_property_map<face_descriptor, Vector>("f:normals", CGAL::NULL_VECTOR).first;

  CGAL::Polygon_mesh_processing::compute_face_normals(m, fnormals);

  for (Edge_index e : m.edges())
  {
    Halfedge_index h = e.halfedge();
    Face_index f1 = m.face(h);
    Face_index f2 = m.face(m.opposite(h));
    const Vector& n1 = fnormals[f1];
    const Vector& n2 = fnormals[f2];
    if (n1*n2 < cos_tolerance)
    {
      longest_protected = std::max(longest_protected,
                                   Vector(m.point(m.source(h)), m.point(m.target(h))).squared_length());
      edge_constrained_map[e] = true;
      num_protected++;
    }
  }

  m.remove_property_map(fnormals);

  std::cout << "Num protected:     " << num_protected << std::endl;
  std::cout << "Longest protected: " << sqrt(longest_protected) << std::endl;

  //------------------------- Remeshing ----------------------------------
  CGAL::Polygon_mesh_processing::isotropic_remeshing(faces(m),
                                                     .03,
                                                     m,
                                                     Parameters::edge_is_constrained_map(edge_constrained_map)
                                                     .protect_constraints(true)
    );

  {
    std::ofstream ofile("remeshed.off");
    ofile << m;
  }
}

Attachment: remeshed.off
Description: model/geomview-off

Attachment: cylinder.off
Description: model/geomview-off