Getting segmentation fault from CGAL::extract_mean_curvature_flow_skeleton

myousefi2016 commented 5 years ago

Issue Details

I'm using CGAL::extract_mean_curvature_flow_skeleton example to extract centerline of a blood vessel 3D surface. But when I run the program, which works fine on default elephant.off example, it shows segmentation fault.

I shared the not capped and capped surfaces here. You can reproduce the segmentation fault error with Case1.off and the code will work flawlessly with Case1-capped.off.

Case1.off Case1-capped.off

Source Code

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/IO/Polyhedron_iostream.h>
#include <CGAL/boost/graph/graph_traits_Polyhedron_3.h>
#include <CGAL/boost/graph/properties_Polyhedron_3.h>
#include <CGAL/extract_mean_curvature_flow_skeleton.h>
#include <CGAL/boost/graph/split_graph_into_polylines.h>
#include <fstream>
#include <boost/foreach.hpp>

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

//only needed for the display of the skeleton as maximal polylines
struct Display_polylines{
  const Skeleton& skeleton;
  std::ofstream& out;
  int polyline_size;
  std::stringstream sstr;
  Display_polylines(const Skeleton& skeleton, std::ofstream& out)
    : skeleton(skeleton), out(out)
  {}
  void start_new_polyline(){
    polyline_size=0;
    sstr.str("");
    sstr.clear();
  }
  void add_node(Skeleton_vertex v){
    ++polyline_size;
    sstr << " " << skeleton[v].point;
  }
  void end_polyline()
  {
    out << polyline_size << sstr.str() << "\n";
  }
};

// This example extracts a medially centered skeleton from a given mesh.
int main(int argc, char* argv[])
{
  // argv[1] Input mesh in .off format
  // argv[2] Output skeleton in .cgal format
  // argv[3] Output point correspondances in .cgal format
  std::ifstream input((argc>1)?argv[1]:"data/elephant.off");
  Polyhedron tmesh;
  input >> tmesh;
  Skeleton skeleton;
  CGAL::extract_mean_curvature_flow_skeleton(tmesh, skeleton);
  std::cout << "Number of vertices of the skeleton: " << boost::num_vertices(skeleton) << "\n";
  std::cout << "Number of edges of the skeleton: " << boost::num_edges(skeleton) << "\n";
  // Output all the edges of the skeleton.
  std::ofstream output((argc>2)?argv[2]:"skel.cgal");
  Display_polylines display(skeleton,output);
  CGAL::split_graph_into_polylines(skeleton, display);
  output.close();
  // Output skeleton points and the corresponding surface points
  output.open((argc>3)?argv[3]:"correspondance.cgal");
  BOOST_FOREACH(Skeleton_vertex v, vertices(skeleton))
    BOOST_FOREACH(vertex_descriptor vd, skeleton[v].vertices)
      output << "2 " << skeleton[v].point << " "
                     << get(CGAL::vertex_point, tmesh, vd)  << "\n";
  return 0;
}

Environment

Operating system (Windows/Mac/Linux, 32/64 bits): Ubuntu 14.04
Compiler: GCC 4.8
Release or debug mode: Release
Specific flags used (if any):
CGAL version: 4.12
Boost version: 1.55
Other libraries versions if used (Eigen, TBB, etc.): Eigen3

Also I run the program by debug flags through gdb and I got this result:

yousefi@oberon:~/skeleton/skeleton/build$ gdb ./skeleton 
GNU gdb (Ubuntu 7.7.1-0ubuntu5~14.04.3) 7.7.1
Copyright (C) 2014 Free Software Foundation, Inc.
License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>
This is free software: you are free to change and redistribute it.
There is NO WARRANTY, to the extent permitted by law.  Type "show copying"
and "show warranty" for details.
This GDB was configured as "x86_64-linux-gnu".
Type "show configuration" for configuration details.
For bug reporting instructions, please see:
<http://www.gnu.org/software/gdb/bugs/>.
Find the GDB manual and other documentation resources online at:
<http://www.gnu.org/software/gdb/documentation/>.
For help, type "help".
Type "apropos word" to search for commands related to "word"...
Reading symbols from ./skeleton...done.
(gdb) run ./Case1.off 
Starting program: /users/yousefi/skeleton/skeleton/build/skeleton ./Case1.off
[Thread debugging using libthread_db enabled]
Using host libthread_db library "/lib/x86_64-linux-gnu/libthread_db.so.1".

Program received signal SIGSEGV, Segmentation fault.
0x00000000007c5381 in CGAL::Mean_curvature_flow_skeletonization<CGAL::Polyhedron_3<CGAL::Simple_cartesian<double>, CGAL::Polyhedron_items_3, CGAL::HalfedgeDS_default, std::allocator<int> >, CGAL::Default, CGAL::Default, CGAL::Default>::compute_voronoi_pole (this=0x7fffffffbf50)
    at /users/yousefi/local/include/CGAL/Mean_curvature_flow_skeletonization.h:1322
1322          p.second->pole = cell_dual[max_neg_i];
(gdb)

I appreciate any help to resolve this issue.

Update

I realized that my 3D surface is not a closed manifold and as a result I capped it and finally it works fine. I think it's good to add this line to make sure the loaded 3D surface does not contain any hole:

  unsigned int nb_holes = 0;
  BOOST_FOREACH(Halfedge_handle h, halfedges(tmesh))
  {
    if(h->is_border())
    {
      std::vector<Facet_handle>  patch_facets;
      std::vector<Vertex_handle> patch_vertices;
        CGAL::Polygon_mesh_processing::triangulate_and_refine_hole(
                  tmesh,
                  h,
                  std::back_inserter(patch_facets),
                  std::back_inserter(patch_vertices),
     CGAL::Polygon_mesh_processing::parameters::vertex_point_map(get(CGAL::vertex_point, tmesh)).
                  geom_traits(Kernel()));
      std::cout << " Number of facets in constructed patch: " << patch_facets.size() << std::endl;
      std::cout << " Number of vertices in constructed patch: " << patch_vertices.size() << std::endl;
      ++nb_holes;
    }
  }

sloriot commented 5 years ago

Could you share the model that made the code crash?

myousefi2016 commented 5 years ago

Could you share the model that made the code crash?

Done.

sloriot commented 5 years ago

From here, you'll read: tmesh is a triangle mesh without borders and having exactly one connected component.

I agree the outcome is a bit harsh but it is the expected behavior. I think in debug mode you'll get a precondition error instead of the segfault. Note that the free functionis_closed() can be used to detect when the input had a boundary or not.

CGAL / cgal