CGAL / cgal

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

CGAL 5.1 extract_mean_curvature_flow_skeleton fails with Bus error: 10 #5084

Closed r03ert0 closed 4 years ago

r03ert0 commented 4 years ago

both_holesVol.off.zip

Issue Details

I'm trying to extract the skeleton of a series of meshes (example attached) [these are brain meshes where all negative curvature regions have been removed, and the result converted to a volume by extrusion along the normal]. extract_mean_curvature_flow_skeleton fails with Bus error: 10, however, the mesh has a single component, only triangles, and no borders.

Source Code

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/extract_mean_curvature_flow_skeleton.h>
#include <CGAL/boost/graph/split_graph_into_polylines.h>
#include <fstream>

#include <CGAL/config.h>
#include <CGAL/Polygon_mesh_processing/repair_polygon_soup.h>
#include <CGAL/Polygon_mesh_processing/orient_polygon_soup.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     K;
typedef K::Point_3                                            Point;
typedef CGAL::Polyhedron_3<K>                                 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;

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::cout<<"CGAL version: "<<CGAL_VERSION_STR<<std::endl;
  std::cout<<"Total arguments: "<<argc<<std::endl;
  std::cout<<"Input mesh: "<<argv[1]<<std::endl;
  std::cout<<"Output skeleton: "<<argv[2]<<std::endl;
  std::cout<<"Output correspondances: "<<argv[3]<<std::endl;

  std::ifstream input((argc>1)?argv[1]:"data/elephant.off");

  std::vector<K::Point_3> points;
  std::vector<std::vector<std::size_t> > triangles;

  std::cout<<"> reading .off input file"<<std::endl;
  CGAL::read_OFF(input, points, triangles);

  std::cout<<"> reparing polygon soup"<<std::endl;
  CGAL::Polygon_mesh_processing::repair_polygon_soup(points, triangles);

  std::cout<<"> orienting polygon soup"<<std::endl;
  CGAL::Polygon_mesh_processing::orient_polygon_soup(points, triangles);

  std::cout<<"> converting polygon soup to polygon mesh"<<std::endl;
  Polyhedron poly;
  CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(points, triangles, poly);

  if ( poly.empty() ) {
    std::cerr << "ERROR: Poly is empty" << std::endl;
    return 1;
  }
  if (!CGAL::is_triangle_mesh(poly)) {
    std::cerr << "ERROR: Poly is not a triangle mesh" << std::endl;
    return 1;
  }
  if (!CGAL::is_closed(poly)) {
    std::cerr << "ERROR: Poly is not closed" << std::endl;
    return 1;
  }

  Skeleton skeleton;
  CGAL::extract_mean_curvature_flow_skeleton(poly, skeleton);

  return 0;
}

Environment

r03ert0 commented 4 years ago

I tried running the individual steps of extract_mean_curvature_flow_skeleton one after the other: contract_geometry, collapse_edges, split_faces and detect_degeneracies. It's the contract_geometry step that throws the error. I also tried starting by running collapse_edges, split_faces and detect_degeneracies in a loop until there's no more collapses, splits or degeneracies. But even then contract_geometry will throw an error.

  Skeleton skeleton;
  Skeletonization mcs(poly);
  std::size_t ncol, nsplit, ndeg;
  do {
    ncol = mcs.collapse_edges();
    nsplit = mcs.split_faces();
    ndeg = mcs.detect_degeneracies();
    std::cout<<"#col:"<<ncol<<" #split:"<<nsplit<<" #deg:"<<ndeg<<std::endl;
  } while(ncol + nsplit + ndeg > 0);
  mcs.contract_geometry();

All my meshes are produced by the same algorithm (remove regions, extrusion). I can't figure out what's that makes that skeletonization doesn't work in some of them...

sloriot commented 4 years ago

Screenshot from 2020-10-19 09-52-50 The two red points (20.8867 25.077 34.0597 and 20.8379 25.2667 33.9117) have been duplicated to handle a non-manifold edge in your mesh. This cause an issue when computing poles in the current code. I'll see if I can easily find a workaround but in any case, the skeleton will not cross the non-manifold edge.

r03ert0 commented 4 years ago

thank you @sloriot!

sloriot commented 4 years ago

You can try the patch in #5088

r03ert0 commented 4 years ago

i works beautifully! 🎉 Thank you very much for the patch. I hope it'll be merged soon 😄

sloriot commented 4 years ago

It was merged today in master and release branches starting from 5.0

afabri commented 4 years ago

Hi @r03ert0 Does the mean curvature flow have something to do with http://neuroanatomy.github.io/ ? If so what should we link on the page Projects Using CGAL.

r03ert0 commented 4 years ago

Hi @afabri! https://neuroanatomy.github.io is the website of my lab. The project where we use the mean curvature flow is this one: https://github.com/r03ert0/foldgraph (as soon as we have a more stable version, we'll move it from my personal repo to the lab's repo). Here's an image showing the workflow we implement with Foldgraph:

foldgraph

And now, we also use the distance information to be able to measure the "width" of each brain fold. We're very excited about the possibilities!

Screenshot 2020-10-20 at 17 09 29
SY524044272-python commented 3 years ago

hi@r03ert0, I want to ask, how do you display the extracted skeleton?I can use this function to generate .cgalfile, but I don't know how to display

r03ert0 commented 3 years ago

hi @SY524044272-python ! I get 2 cgal files from the skeleton sample code: "skel-poly.cgal" and "correspondance-poly.cgal". For getting the skeleton/wireframe U use the skel-poly file. Each row is one curve, the first value is the number of vertices in the curve, and then follow the 3d coordinates of the vertices. After that, the "correspondance-poly" file links each point in the skeleton curve with its corresponding point in the original mesh. For plotting in 3d I used custom code, but it should be possible to convert the curves to PLY format or something, and render in Blender (you can get some very cool renders and animations!). For the flat plots, I used some custom code for the flattening (which is quite specific to our brain imaging pipeline), and then just plotted using Matplotlib : ). In a previous experiment, I used t-SNE for flattening the graphs, and it works more or less :P

SY524044272-python commented 3 years ago

hello Roberto Toro. I'm sorry to see your reply very late. I already know how to use these two files。 But I still have a few questions to ask you. I use SolidWorks to design a chain tooth structure and output it as STL file, then write Python code to convert the STL file into the off file. When I use extract mean curvature flow Skeleton, there was an error in the program. It's like the picture below.I searched GitHub for this problem and found that there were also people who ran into it.

This is the solution, but I didn't understand what he said

My chain tooth file is in the attachment. Can you help me see what the problem is? thank you very much!!!!!!

I'm sorry to see your reply very late

At 2020-12-10 17:58:13, "Roberto Toro" @.***> wrote:

hi @SY524044272-python ! I get 2 cgal files from the skeleton sample code: "skel-poly.cgal" and "correspondance-poly.cgal". For getting the skeleton/wireframe U use the skel-poly file. Each row is one curve, the first value is the number of vertices in the curve, and then follow the 3d coordinates of the vertices. After that, the "correspondance-poly" file links each point in the skeleton curve with its corresponding point in the original mesh. For plotting in 3d I used custom code, but it should be possible to convert the curves to PLY format or something, and render in Blender (you can get some very cool renders and animations!). For the flat plots, I used some custom code for the flattening (which is quite specific to our brain imaging pipeline), and then just plotted using Matplotlib : ). In a previous experiment, I used t-SNE for flattening the graphs, and it works more or less :P

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.