Mesh cleaning (duplicated vertices removal)

[stla]

Hello,

When I do an Advanced Front Surface reconstruction, I end up with duplicated vertices and the triangles indices are simply 1, 2, 3 - 4, 5, 6 - ...

The result is sent to R and I do the mesh cleaning in R, with the Rvcg package. It removes the duplicated vertices and recalculates the faces indices. Is it possible to do that with CGAL?

Below is my code but it uses the Rcpp library, this is the library which allows to combine R and C++.

Rcpp::List AFSreconstruction_cpp(Rcpp::NumericMatrix pts) {
  const size_t npoints = pts.nrow();
  std::vector<Point3> points(npoints);
  for(size_t i = 0; i < npoints; i++) {
    points[i] = Point3(pts(i, 0), pts(i, 1), pts(i, 2));
  }

  AFS_triangulation3 dt(points.begin(), points.end());
  AFS_reconstruction reconstruction(dt);
  reconstruction.run();
  const AFS_Tds2& tds = reconstruction.triangulation_data_structure_2();

  Eigen::MatrixXd vertices(4, 0);
  unsigned counter = 0;
  for(AFS_Tds2::Face_iterator fit = tds.faces_begin(); fit != tds.faces_end();
      ++fit) {
    if(reconstruction.has_on_surface(fit)) {
      counter++;
      AFS_triangulation3::Facet f = fit->facet();
      AFS_triangulation3::Cell_handle ch = f.first;
      int ci = f.second;
      Point3 points[3];
      for(int i = 0, j = 0; i < 4; i++) {
        if(ci != i) {
          points[j] = ch->vertex(i)->point();
          j++;
        }
      }
      for(size_t k = 0; k < 3; k++) {
        const Point3 p = points[k];
        Eigen::VectorXd w(4);
        w << p.x(), p.y(), p.z(), 1.0;
        vertices.conservativeResize(Eigen::NoChange, vertices.cols() + 1);
        vertices.rightCols(1) = w;
      }
    }
  }
  Rcpp::IntegerVector vtriangles(3 * counter);
  for(size_t i = 0; i < 3 * counter; i++) {
    vtriangles(i) = i + 1;
  }
  vtriangles.attr("dim") = Rcpp::Dimension(3, counter);
  Rcpp::IntegerMatrix triangles = Rcpp::as<Rcpp::IntegerMatrix>(vtriangles);

  return Rcpp::List::create(Rcpp::Named("vertices") = vertices,
                            Rcpp::Named("triangles") = triangles);
}

[sloriot]

You need to assign to each vertex an index that you will then use when describing a face.

You can use the following code to set and get the id of a vertex per face (did not try to compile it).

int i=0;
std::unordered_map<AFS_Tds2::Vertex_Handle, int> v2i;
auto get_face_id = [&v2i, &i](AFS_Tds2::Vertex_Handle v)
{
  auto insert_res = v2i.insert( std::make_pair(v, i));
  if (insert_res.second) ++i;
  return insert_res.first->second;
};

....

int vid = get_face_id(ch->vertex(i));

[sloriot]

Even simpler, adapt the following example: https://doc.cgal.org/latest/Advancing_front_surface_reconstruction/Advancing_front_surface_reconstruction_2reconstruction_surface_mesh_8cpp-example.html

Vertices are provided by batch and you have a call to operator() per face.

[stla]

Thanks. When I deal with C++ I'm playing beyond my skills. Here I see a vector of Facet which is declared (std::vector<Facet> facets;) but it is used nowhere. Is it normal?

[sloriot]

Indeed, I guess at first @afabri wanted to fill a soup but decided to fill a Surface_mesh instead. @afabri would you mind fixing it?

If you want to use the two vectors instead of the mesh, simply push_back points and faces in the constructor and operator=.