CGAL / cgal

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

Different results are obtained in different runs with the same input of CGAL::Mean_curvature_flow_skeletonization. #7511

Closed qq775193759 closed 1 year ago

qq775193759 commented 1 year ago

Issue Details

We demonstrate a bug (about stability) of CGAL::Mean_curvature_flow_skeletonization. We use Mean_curvature_flow_skeletonization to generate skeleton in our work and obtain different results in different runs, which is adverse to testing and optimizing our work. Therefore, we think it is necessary to fix the bug.

We notice that the issue #5726 has proposed a similar issue. A comment in issue #5726 is "It is a known issue that it is not deterministic, probably due to Polyhedron handles being hashes or put in a map". We find the cause of the uncertainty is the less operation, which uses a not deterministic key to compare two edges.

There are 4 key steps in Skeletonization: contract_geometry(), collapse_edges(), split_faces() and detect_degeneracies(). We find collapse_edges() is the only cause (by testing with parts of code).

The whole class CGAL::Mean_curvature_flow_skeletonization is in the file CGAL/Mean_curvature_flow_skeletonization.h: In Mean_curvature_flow_skeletonization::collapse_short_edges(): it use std::set. The elements are sorted by the less operation. I don't know the implement of the less operation (someone could help me). I print the order of elements and they are different in different runs. Therefore, it must be the cause.

Based on the observation, we fix the code by a simple way. We overrid the less operation of edge_descriptor and obtain deterministic results. The change of CGAL code is the following:

template <class T>
struct IdLess {
  bool operator()(const T& x, const T& y) const { return x.id() < y.id(); }
};

template <class TriangleMesh, class Traits_, class VertexPointMap_, class SolverTraits_>
std::size_t Mean_curvature_flow_skeletonization<TriangleMesh, Traits_, VertexPointMap_,
                                                SolverTraits_>::collapse_short_edges() {
  std::size_t cnt = 0, prev_cnt = 0;

  std::set<edge_descriptor, IdLess<edge_descriptor>> edges_to_collapse,
      non_topologically_valid_collapses;

Source Code (The input file and the output are in skeleton_debug.zip)

#include <CGAL/Mean_curvature_flow_skeletonization.h>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Surface_mesh.h>

#include <fstream>

typedef CGAL::Simple_cartesian<double> Kernel;
typedef Kernel::Point_3 Point;
typedef CGAL::Surface_mesh<Point> Triangle_mesh;
typedef boost::graph_traits<Triangle_mesh>::vertex_descriptor vertex_descriptor;
typedef CGAL::Mean_curvature_flow_skeletonization<Triangle_mesh>
    Skeletonization;
typedef Skeletonization::Skeleton Skeleton;
typedef Skeleton::vertex_descriptor Skeleton_vertex;
typedef Skeleton::edge_descriptor Skeleton_edge;

void MeshToSkeleton(const std::string filename) {
  Triangle_mesh tmesh;
  CGAL::IO::read_OBJ(filename, tmesh);
  if (!CGAL::is_triangle_mesh(tmesh)) {
    std::cout << "Input geometry is not triangulated." << std::endl;
    return;
  }
  static int file_number = 0;
  Skeleton skeleton;
  Skeletonization mcs(tmesh);
  mcs.contract_until_convergence();
  mcs.convert_to_skeleton(skeleton);

  std::cout << "running process number: " << file_number << std::endl;
  std::cout << "Number of vertices of the skeleton: "
            << boost::num_vertices(skeleton) << std::endl;
  std::cout << "Number of edges of the skeleton: " << boost::num_edges(skeleton)
            << std::endl
            << std::endl;
  // Output all the edges of the skeleton.
  std::ofstream output(std::to_string(file_number) + "_skel.txt");
  for (Skeleton_edge e : CGAL::make_range(edges(skeleton))) {
    const Point& s = skeleton[source(e, skeleton)].point;
    const Point& t = skeleton[target(e, skeleton)].point;
    output << "2 " << s << " " << t << "\n";
  }
  output.close();
  file_number++;
}

int main(int argc, char* argv[]) {
  const std::string filename = "debug.obj";
  MeshToSkeleton(filename);
  MeshToSkeleton(filename);
  MeshToSkeleton(filename);
  return EXIT_SUCCESS;
}

Command Line Output

running process number: 0 Number of vertices of the skeleton: 215 Number of edges of the skeleton: 214

running process number: 1 Number of vertices of the skeleton: 223 Number of edges of the skeleton: 222

running process number: 2 Number of vertices of the skeleton: 220 Number of edges of the skeleton: 219

Environment

afabri commented 1 year ago

Thank you for pointing out the problem as well as the fix.

lrineau commented 1 year ago

Fixed by https://github.com/CGAL/cgal/pull/7551.