Bump libigl to f68f01c2 and add bindings for new/missing functions

CMakeLists.txt

@@ -60,7 +60,7 @@ option(LIBIGL_CYCODEBASE           "Build igl::cycodebase bindings"        ON)
 FetchContent_Declare(
   libigl
   GIT_REPOSITORY https://github.com/libigl/libigl.git
-  GIT_TAG 0cbcaf886314ff331402c38fbb04cc4de8d7c033
+  GIT_TAG f68f01c23ac82351e2543ce8c529b6d4c89fdf16
 )
 FetchContent_MakeAvailable(libigl)
 

src/face_areas.cpp

@@ -17,6 +17,14 @@ namespace pyigl
     igl::face_areas(V, T, A);
     return A;
   }
+
+  auto face_areas_intrinsic(
+    const nb::DRef<const Eigen::MatrixXN>& L)
+  {
+    Eigen::MatrixXN A;
+    igl::face_areas(L, A);
+    return A;
+  }
 }
 
 void bind_face_areas(nb::module_ &m)
@@ -27,6 +35,13 @@ void bind_face_areas(nb::module_ &m)
     R"(Constructs a list of face areas of faces opposite each index in a tet list
 @param[in] V  #V by 3 list of mesh vertex positions
 @param[in] T  #T by 3 list of tet mesh indices into V
+@param[out] A   #T by 4 list of face areas corresponding to faces opposite vertices
+    0,1,2,3)");
+  m.def("face_areas", &pyigl::face_areas_intrinsic,
+    "L"_a,
+    R"(Compute tet-mesh face areas from edge lengths.
+@param[in] L  #T by 6 list of tet-mesh edge lengths corresponding to edges
+    [1,2],[2,0],[0,1],[3,0],[3,1],[3,2]
 @param[out] A   #T by 4 list of face areas corresponding to faces opposite vertices
     0,1,2,3)");
 }

src/per_face_normals.cpp

@@ -20,6 +20,15 @@ namespace pyigl
     igl::per_face_normals(V,F,Z,N);
     return N;
   }
+  // Wrapper for per_face_normals_stable function
+  auto per_face_normals_stable(
+    const nb::DRef<const Eigen::MatrixXN> &V,
+    const nb::DRef<const Eigen::MatrixXI> &F)
+  {
+    Eigen::MatrixXN N;
+    igl::per_face_normals_stable(V,F,N);
+    return N;
+  }
   // Wrapper for per_face_normals function
   auto per_face_normals_VIC(
     const nb::DRef<const Eigen::MatrixXN> &V,
@@ -69,4 +78,17 @@ R"(Compute face normals via vertex position list, polygon stream
 @param[out] FF  #I by 3 list of triangle indices into rows of VV
 @param[out] J  #I list of indices into original polygons)"
     );
+  m.def(
+    "per_face_normals_stable",
+    &pyigl::per_face_normals_stable,
+    "V"_a,
+    "F"_a,
+R"(Special version of per_face_normals where the order of the face indices is
+guaranteed not to affect the output (i.e. cross products are accumulated in a
+consistent order). Degenerate faces are given a zero normal.
+
+@param[in] V  #V by 3 eigen Matrix of mesh vertex 3D positions
+@param[in] F  #F by 3 eigen Matrix of face (triangle) indices
+@param[out] N  #F by 3 eigen Matrix of mesh face (triangle) 3D normals)"
+    );
 }

src/simplex_simplex_squared_distance.cpp

@@ -0,0 +1,44 @@
+#include "default_types.h"
+#include <igl/simplex_simplex_squared_distance.h>
+#include <nanobind/nanobind.h>
+#include <nanobind/ndarray.h>
+#include <nanobind/eigen/dense.h>
+#include <nanobind/stl/tuple.h>
+
+namespace nb = nanobind;
+using namespace nb::literals;
+
+namespace pyigl
+{
+  // Wrapper for simplex_simplex_squared_distance function
+  auto simplex_simplex_squared_distance(
+    const nb::DRef<const Eigen::MatrixXN> &V1,
+    const nb::DRef<const Eigen::MatrixXN> &V2)
+  {
+    Numeric sqrdDist;
+    Eigen::VectorXN B1;
+    Eigen::VectorXN B2;
+    igl::simplex_simplex_squared_distance(V1, V2, sqrdDist, B1, B2);
+    return std::make_tuple(sqrdDist, B1, B2);
+  }
+}
+
+// Bind the wrapper to the Python module
+void bind_simplex_simplex_squared_distance(nb::module_ &m)
+{
+  m.def(
+    "simplex_simplex_squared_distance",
+    &pyigl::simplex_simplex_squared_distance,
+    "V1"_a,
+    "V2"_a,
+R"(Find the squared distance between closest points on simplices with corners
+V1 and V2, respectively. V1 and V2 don't have to be the same simplex size,
+but they must have the same number of columns (dimension). This function
+works recursively.
+
+@param[in] V1  #V1 by dim list of simplex corners
+@param[in] V2  #V2 by dim list of simplex corners
+@param[out] sqrdDist  squared distance between closest points on simplices
+@param[out] B1  #V1 list of barycentric coordinates of closest point on simplex 1
+@param[out] B2  #V2 list of barycentric coordinates of closest point on simplex 2)");
+}

src/triangle/remesh_at_points.cpp

@@ -0,0 +1,55 @@
+#include "default_types.h"
+#include <igl/triangle/remesh_at_points.h>
+#include <nanobind/nanobind.h>
+#include <nanobind/ndarray.h>
+#include <nanobind/eigen/dense.h>
+#include <nanobind/stl/tuple.h>
+
+namespace nb = nanobind;
+using namespace nb::literals;
+
+namespace pyigl
+{
+  // Wrapper for triangle::remesh_at_points function
+  auto remesh_at_points(
+    const nb::DRef<const Eigen::MatrixXN> &V,
+    const nb::DRef<const Eigen::MatrixXI> &F,
+    const nb::DRef<const Eigen::MatrixXN> &B,
+    const nb::DRef<const Eigen::VectorXI> &FI)
+  {
+    Eigen::MatrixXN VV;
+    Eigen::MatrixXI FF;
+    Eigen::VectorXI J;
+    Eigen::VectorXI K;
+    igl::triangle::remesh_at_points(V, F, B, FI, VV, FF, J, K);
+    return std::make_tuple(VV, FF, J, K);
+  }
+}
+
+void bind_remesh_at_points(nb::module_ &m)
+{
+  m.def(
+    "remesh_at_points",
+    &pyigl::remesh_at_points,
+    "V"_a,
+    "F"_a,
+    "B"_a,
+    "FI"_a,
+R"(Given a set of unique points on a mesh specified by barycentric coordinates
+and a triangle index, remesh the mesh to include these points as vertices.
+Barycentric coordinates should be non-negative and sum to 1. Vertex-points are
+not inserted (but tracked in K). Edge-points are inserted and the edge is split
+on all incident faces. Face-points are inserted. The output preserves the
+(non-)manifoldness of the input. The input mesh can be in 3D (or any dimension).
+
+@param[in] V  #V by dim list of mesh vertex positions
+@param[in] F  #F by 3 list of triangle indices into rows of V
+@param[in] B  #B by 3 list of barycentric coordinates, ith row are coordinates
+              of ith sampled point in face FI(i)
+@param[in] FI  #B list of indices into F
+@param[out] VV  #VV by dim list of mesh vertex positions (top #V rows is always
+                V; bottom rows are new vertices)
+@param[out] FF  #FF by 3 list of triangle indices into rows of VV
+@param[out] J  #FF list of indices into F
+@param[out] K  #B list of indices into VV)");
+}

tests/test_all.py

@@ -142,6 +142,9 @@ def test_normals_and_distances():
     FN,_,_,_ = igl.per_face_normals(V,I,C)
     FN = igl.per_face_normals(V,F)
     FN = igl.per_face_normals(V,F,Z=np.array([0,0,1],dtype=np.float64))
+    FNs = igl.per_face_normals_stable(V,F)
+    assert FNs.shape == FN.shape
+    np.testing.assert_allclose(np.abs((FNs*FN).sum(axis=1)), 1.0, atol=1e-9)
     VN = igl.per_vertex_normals(V,F)
     VN = igl.per_vertex_normals(V,F,      weighting=igl.PER_VERTEX_NORMALS_WEIGHTING_TYPE_UNIFORM)
     VN = igl.per_vertex_normals(V,F,      weighting=igl.PER_VERTEX_NORMALS_WEIGHTING_TYPE_AREA)
@@ -625,6 +628,9 @@ def test_misc():
     theta, cos_theta = igl.dihedral_angles(V,T)
     L = igl.edge_lengths(V,T)
     A = igl.face_areas(V,T)
+    # intrinsic overload: face areas from tet edge lengths alone
+    A_intrinsic = igl.face_areas(L)
+    np.testing.assert_allclose(A_intrinsic, A, atol=1e-9)
     theta, cos_theta = igl.dihedral_angles_intrinsic(L,A)
     D = igl.all_pairs_distances(V,V,squared=False)
     D = igl.all_pairs_distances(V,V,squared=True)
@@ -1423,6 +1429,61 @@ def test_new_triangle_algorithms():
     assert V_ref.shape[0] >= V_tri.shape[0]
 
 
+def test_remesh_at_points():
+    # Single triangle in 3D
+    V = np.array([[0.0, 0.0, 0.0],
+                  [1.0, 0.0, 0.0],
+                  [0.0, 1.0, 0.0]], dtype=np.float64)
+    F = np.array([[0, 1, 2]], dtype=np.int64)
+    # One interior (face) point at the centroid and one on an edge midpoint
+    B = np.array([[1/3, 1/3, 1/3],
+                  [0.5, 0.5, 0.0]], dtype=np.float64)
+    FI = np.array([0, 0], dtype=np.int64)
+    VV, FF, J, K = igl.triangle.remesh_at_points(V, F, B, FI)
+    # Output vertices keep the input as the top rows
+    assert VV.shape[1] == 3
+    assert VV.shape[0] >= V.shape[0]
+    np.testing.assert_allclose(VV[:V.shape[0]], V)
+    # New faces index into VV and map back to original faces via J
+    assert FF.shape[1] == 3
+    assert J.shape[0] == FF.shape[0]
+    assert FF.max() < VV.shape[0]
+    assert J.max() < F.shape[0]
+    # K tracks each sampled point's index into VV
+    assert K.shape[0] == B.shape[0]
+    assert K.max() < VV.shape[0]
+    # The centroid face-point must have been inserted as a real vertex
+    np.testing.assert_allclose(VV[K[0]], B[0] @ V)
+
+
+def test_simplex_simplex_squared_distance():
+    # Two parallel segments in 3D separated by distance 1 along z
+    V1 = np.array([[0.0, 0.0, 0.0],
+                   [1.0, 0.0, 0.0]], dtype=np.float64)
+    V2 = np.array([[0.0, 0.0, 1.0],
+                   [1.0, 0.0, 1.0]], dtype=np.float64)
+    sqrD, B1, B2 = igl.simplex_simplex_squared_distance(V1, V2)
+    assert np.isclose(sqrD, 1.0)
+    # Barycentric coordinates: one per corner of each simplex, summing to 1
+    assert B1.shape[0] == V1.shape[0]
+    assert B2.shape[0] == V2.shape[0]
+    assert np.isclose(B1.sum(), 1.0)
+    assert np.isclose(B2.sum(), 1.0)
+    # Reconstructed closest points should be 1 apart
+    P1 = B1 @ V1
+    P2 = B2 @ V2
+    assert np.isclose(np.linalg.norm(P1 - P2), 1.0)
+
+    # Point-to-triangle: point directly above a triangle's centroid
+    T = np.array([[0.0, 0.0, 0.0],
+                  [1.0, 0.0, 0.0],
+                  [0.0, 1.0, 0.0]], dtype=np.float64)
+    P = np.array([[1/3, 1/3, 2.0]], dtype=np.float64)
+    sqrD2, Bp, Bt = igl.simplex_simplex_squared_distance(P, T)
+    assert np.isclose(sqrD2, 4.0)
+    assert np.isclose(Bt.sum(), 1.0)
+
+
 def test_lexicographic_triangulation():
     # Simple square: 4 points in general position
     P = np.array([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]], dtype=np.float64)