203_CurvatureDirections.py 1.3 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556
  1. import igl
  2. std::string filename = "../shared/fertility.off";
  3. if(argc>1)
  4. {
  5. filename = argv[1];
  6. }
  7. // Load a mesh in OFF format
  8. igl::read_triangle_mesh(filename, V, F);
  9. // Alternative discrete mean curvature
  10. MatrixXd HN;
  11. SparseMatrix<double> L,M,Minv;
  12. igl::cotmatrix(V,F,L);
  13. igl::massmatrix(V,F,igl::MASSMATRIX_TYPE_VORONOI,M);
  14. igl::invert_diag(M,Minv);
  15. // Laplace-Beltrami of position
  16. HN = -Minv*(L*V);
  17. // Extract magnitude as mean curvature
  18. VectorXd H = HN.rowwise().norm();
  19. // Compute curvature directions via quadric fitting
  20. MatrixXd PD1,PD2;
  21. VectorXd PV1,PV2;
  22. igl::principal_curvature(V,F,PD1,PD2,PV1,PV2);
  23. // mean curvature
  24. H = 0.5*(PV1+PV2);
  25. igl::viewer::Viewer viewer;
  26. viewer.data.set_mesh(V, F);
  27. // Compute pseudocolor
  28. MatrixXd C;
  29. igl::parula(H,true,C);
  30. viewer.data.set_colors(C);
  31. // Average edge length for sizing
  32. const double avg = igl::avg_edge_length(V,F);
  33. // Draw a blue segment parallel to the minimal curvature direction
  34. const RowVector3d red(0.8,0.2,0.2),blue(0.2,0.2,0.8);
  35. viewer.data.add_edges(V + PD1*avg, V - PD1*avg, blue);
  36. // Draw a red segment parallel to the maximal curvature direction
  37. viewer.data.add_edges(V + PD2*avg, V - PD2*avg, red);
  38. // Hide wireframe
  39. viewer.core.show_lines = false;
  40. viewer.launch();
  41. }