203_CurvatureDirections.py 2.0 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879
  1. #!/usr/bin/env python
  2. #
  3. # This file is part of libigl, a simple c++ geometry processing library.
  4. #
  5. # Copyright (C) 2017 Sebastian Koch <[email protected]> and Daniele Panozzo <[email protected]>
  6. #
  7. # This Source Code Form is subject to the terms of the Mozilla Public License
  8. # v. 2.0. If a copy of the MPL was not distributed with this file, You can
  9. # obtain one at http://mozilla.org/MPL/2.0/.
  10. import sys, os
  11. # Add the igl library to the modules search path
  12. sys.path.insert(0, os.getcwd() + "/../")
  13. import pyigl as igl
  14. from shared import TUTORIAL_SHARED_PATH, check_dependencies
  15. dependencies = ["glfw"]
  16. check_dependencies(dependencies)
  17. V = igl.eigen.MatrixXd()
  18. F = igl.eigen.MatrixXi()
  19. igl.read_triangle_mesh(TUTORIAL_SHARED_PATH + "fertility.off", V, F)
  20. # Alternative discrete mean curvature
  21. HN = igl.eigen.MatrixXd()
  22. L = igl.eigen.SparseMatrixd()
  23. M = igl.eigen.SparseMatrixd()
  24. Minv = igl.eigen.SparseMatrixd()
  25. igl.cotmatrix(V, F, L)
  26. igl.massmatrix(V, F, igl.MASSMATRIX_TYPE_VORONOI, M)
  27. igl.invert_diag(M, Minv)
  28. # Laplace-Beltrami of position
  29. HN = -Minv * (L * V)
  30. # Extract magnitude as mean curvature
  31. H = HN.rowwiseNorm()
  32. # Compute curvature directions via quadric fitting
  33. PD1 = igl.eigen.MatrixXd()
  34. PD2 = igl.eigen.MatrixXd()
  35. PV1 = igl.eigen.MatrixXd()
  36. PV2 = igl.eigen.MatrixXd()
  37. igl.principal_curvature(V, F, PD1, PD2, PV1, PV2)
  38. # Mean curvature
  39. H = 0.5 * (PV1 + PV2)
  40. viewer = igl.glfw.Viewer()
  41. viewer.data().set_mesh(V, F)
  42. # Compute pseudocolor
  43. C = igl.eigen.MatrixXd()
  44. igl.parula(H, True, C)
  45. viewer.data().set_colors(C)
  46. # Average edge length for sizing
  47. avg = igl.avg_edge_length(V, F)
  48. # Draw a blue segment parallel to the minimal curvature direction
  49. red = igl.eigen.MatrixXd([[0.8, 0.2, 0.2]])
  50. blue = igl.eigen.MatrixXd([[0.2, 0.2, 0.8]])
  51. viewer.data().add_edges(V + PD1 * avg, V - PD1 * avg, blue)
  52. # Draw a red segment parallel to the maximal curvature direction
  53. viewer.data().add_edges(V + PD2 * avg, V - PD2 * avg, red)
  54. # Hide wireframe
  55. viewer.data().show_lines = False
  56. viewer.launch()