ADD curvature

Author: petrasvestartas

docs/examples/curvature.py

@@ -2,49 +2,79 @@
 import compas_libigl as igl
 from compas.colors import Color
 from compas.datastructures import Mesh
-from compas.geometry import Line
-from compas.geometry import Point
-from compas.geometry import Vector
+from compas.geometry import Line, Point, Vector, Scale, Rotation
 from compas_viewer import Viewer
 
 # ==============================================================================
 # Input geometry
 # ==============================================================================
 
 mesh = Mesh.from_obj(compas.get("tubemesh.obj"))
-
 trimesh = mesh.copy()
 trimesh.quads_to_triangles()
 
 # ==============================================================================
-# curvature
+# Curvature
 # ==============================================================================
 
-curvature = igl.trimesh_gaussian_curvature(trimesh.to_vertices_and_faces())
+vertices, faces = trimesh.to_vertices_and_faces()
+gaussian_curvature = igl.trimesh_gaussian_curvature((vertices, faces))
+PD1, PD2, PV1, PV2 = igl.trimesh_principal_curvature((vertices, faces))
+
+# Normalize curvature values for better visualization
+max_gaussian = max(abs(k) for k in gaussian_curvature)
+max_principal = max(max(abs(k) for k in PV1), max(abs(k) for k in PV2))
+
+# Scale factors for visualization
+gaussian_scale = 0.5 / max_gaussian if max_gaussian > 0 else 0.5
+principal_scale = 0.3 / max_principal if max_principal > 0 else 0.3
 
 # ==============================================================================
-# Visualisation
+# Visualization
 # ==============================================================================
 
 viewer = Viewer(width=1600, height=900)
-# viewer.view.camera.position = [1, -6, 2]
-# viewer.view.camera.look_at([1, 1, 1])
 
-viewer.scene.add(mesh, opacity=0.7, show_points=False)
+# Add base mesh
+viewer.scene.add(
+    mesh,
+    opacity=0.7,
+    show_points=False
+)
 
+# Visualize curvature at each vertex
 for vertex in mesh.vertices():
     if mesh.is_vertex_on_boundary(vertex):
         continue
-
+        
+    # Get vertex position and normal
     point = Point(*mesh.vertex_coordinates(vertex))
     normal = Vector(*mesh.vertex_normal(vertex))
-    c = curvature[vertex]
-    normal.scale(10 * c)
-
+    
+    # Gaussian curvature visualization (vertical lines)
+    k = gaussian_curvature[vertex]
+    scaled_normal = normal.scaled(gaussian_scale * k)
     viewer.scene.add(
-        Line(point, point + normal),
-        linecolor=(Color.red() if c > 0 else Color.blue()),
-        linewidth=2,
+        Line(point, point + scaled_normal*100),
+        linecolor=Color.blue(),
+        linewidth=3
     )
+    
+    # Principal curvature visualization (cross lines)
+    pd1 = Vector(*PD1[vertex]).scaled(2*principal_scale * abs(PV1[vertex]))
+    pd2 = Vector(*PD2[vertex]).scaled(2*principal_scale * abs(PV2[vertex]))
+    
+    # Create lines for principal directions
+    viewer.scene.add(
+        Line(point - pd1, point + pd1),
+        linecolor=Color.black(),
+        linewidth=4
+    )
+    viewer.scene.add(
+        Line(point - pd2, point + pd2),
+        linecolor=Color.black(),
+        linewidth=2
+    )
+
 
 viewer.show()

src/compas_libigl/__init__.py

@@ -2,6 +2,7 @@
 import compas
 from ._nanobind import add, __doc__
 from .boundaries import trimesh_boundaries
+from .curvature import trimesh_gaussian_curvature, trimesh_principal_curvature
 
 # from .boundaries import trimesh_boundaries
 # from .curvature import trimesh_gaussian_curvature, trimesh_principal_curvature
@@ -85,7 +86,7 @@ def get_armadillo():
 
 __all_plugins__ = [
     "compas_libigl._nanobindcompas_libigl._boundaries",
-    # "compas_libigl.curvature",
+    "compas_libigl.curvature",
     # "compas_libigl.geodistance",
     # "compas_libigl.intersections",
     # "compas_libigl.isolines",
@@ -103,6 +104,8 @@ def get_armadillo():
     add,
     __doc__,
     trimesh_boundaries,
+    trimesh_gaussian_curvature,
+    trimesh_principal_curvature,
     # "get",
     # "get_beetle",
     # "get_armadillo",

src/compas_libigl/curvature.py

@@ -0,0 +1,68 @@
+import numpy as np
+
+from compas_libigl import _curvature
+
+
+def trimesh_gaussian_curvature(M):
+    """Compute the discrete gaussian curvature of a triangle mesh.
+
+    Parameters
+    ----------
+    M : tuple[list[list[float]], list[list[int]]]
+        A mesh represented by a list of vertices and a list of faces.
+
+    Returns
+    -------
+    list[float]
+        The gaussian curvature per vertex.
+
+    Examples
+    --------
+    >>> import compas
+    >>> import compas_libigl as libigl
+    >>> from compas.datastructures import Mesh
+    >>> mesh = Mesh.from_off(compas.get("tubemesh.off"))
+    >>> mesh.quads_to_triangles()
+    >>> M = mesh.to_vertices_and_faces()
+    >>> K = libigl.trimesh_gaussian_curvature(M)
+    >>> len(K) == len(mesh.vertices)
+    True
+
+    """
+    V, F = M
+    V = np.asarray(V, dtype=np.float64)
+    F = np.asarray(F, dtype=np.int32)
+    return _curvature.trimesh_gaussian_curvature(V, F)
+
+
+def trimesh_principal_curvature(M):
+    """Compute the principal curvatures and directions of a triangle mesh.
+
+    Parameters
+    ----------
+    M : tuple[list[list[float]], list[list[int]]]
+        A mesh represented by a list of vertices and a list of faces.
+
+    Returns
+    -------
+    tuple[list[list[float]], list[list[float]], list[float], list[float]]
+        The principal directions (PD1, PD2) and principal values (PV1, PV2).
+
+    Examples
+    --------
+    >>> import compas
+    >>> import compas_libigl as libigl
+    >>> from compas.datastructures import Mesh
+    >>> mesh = Mesh.from_off(compas.get("tubemesh.off"))
+    >>> mesh.quads_to_triangles()
+    >>> M = mesh.to_vertices_and_faces()
+    >>> PD1, PD2, PV1, PV2 = libigl.trimesh_principal_curvature(M)
+    >>> len(PV1) == len(mesh.vertices)
+    True
+
+    """
+    V, F = M
+    V = np.asarray(V, dtype=np.float64)
+    F = np.asarray(F, dtype=np.int32)
+    PD1, PD2, PV1, PV2 = _curvature.trimesh_principal_curvature(V, F)
+    return PD1.tolist(), PD2.tolist(), PV1.tolist(), PV2.tolist()

src/curvature.cpp

@@ -1,13 +1,42 @@
 #include "curvature.hpp"
 
-void function(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F) {}
+using namespace Eigen;
+
+Eigen::VectorXd trimesh_gaussian_curvature(const compas::RowMatrixXd& V, const compas::RowMatrixXi& F) {
+    Eigen::VectorXd C;
+    igl::gaussian_curvature(V, F, C);
+    return C;
+}
+
+std::tuple<Eigen::MatrixXd, Eigen::MatrixXd, Eigen::VectorXd, Eigen::VectorXd>
+trimesh_principal_curvature(const compas::RowMatrixXd& V, const compas::RowMatrixXi& F) {
+    if (F.cols() != 3) {
+        std::cerr << "Error: Principal curvature requires triangular faces." << std::endl;
+        return std::make_tuple(Eigen::MatrixXd(), Eigen::MatrixXd(),
+                             Eigen::VectorXd(), Eigen::VectorXd());
+    }
+
+    Eigen::MatrixXd PD1, PD2;
+    Eigen::VectorXd PV1, PV2;
+    std::vector<int> bad_vertices;
+
+    igl::principal_curvature(V, F, PD1, PD2, PV1, PV2, bad_vertices);
+
+    return std::make_tuple(PD1, PD2, PV1, PV2);
+}
 
 NB_MODULE(_curvature, m) {
+    m.def(
+        "trimesh_gaussian_curvature",
+        &trimesh_gaussian_curvature,
+        "Compute the discrete gaussian curvature of a triangle mesh.",
+        "V"_a, "F"_a
+    );
 
     m.def(
-        "function_name",
-        &function,
-        "Description.",
-        "V"_a,
-        "F"_a);
+        "trimesh_principal_curvature",
+        &trimesh_principal_curvature,
+        "Compute the principal curvatures and directions of a triangle mesh.",
+        "V"_a, "F"_a
+    );
 }