diff --git a/CMakeLists.txt b/CMakeLists.txt index d40966b7..9ebcd4b7 100644 --- a/CMakeLists.txt +++ b/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) diff --git a/src/face_areas.cpp b/src/face_areas.cpp index c39b4ad0..355548e9 100644 --- a/src/face_areas.cpp +++ b/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& 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)"); } diff --git a/src/per_face_normals.cpp b/src/per_face_normals.cpp index c9908b5a..b82c536e 100644 --- a/src/per_face_normals.cpp +++ b/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 &V, + const nb::DRef &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 &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)" + ); } diff --git a/src/simplex_simplex_squared_distance.cpp b/src/simplex_simplex_squared_distance.cpp new file mode 100644 index 00000000..9fdbd280 --- /dev/null +++ b/src/simplex_simplex_squared_distance.cpp @@ -0,0 +1,44 @@ +#include "default_types.h" +#include +#include +#include +#include +#include + +namespace nb = nanobind; +using namespace nb::literals; + +namespace pyigl +{ + // Wrapper for simplex_simplex_squared_distance function + auto simplex_simplex_squared_distance( + const nb::DRef &V1, + const nb::DRef &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)"); +} diff --git a/src/triangle/remesh_at_points.cpp b/src/triangle/remesh_at_points.cpp new file mode 100644 index 00000000..b660c9fe --- /dev/null +++ b/src/triangle/remesh_at_points.cpp @@ -0,0 +1,55 @@ +#include "default_types.h" +#include +#include +#include +#include +#include + +namespace nb = nanobind; +using namespace nb::literals; + +namespace pyigl +{ + // Wrapper for triangle::remesh_at_points function + auto remesh_at_points( + const nb::DRef &V, + const nb::DRef &F, + const nb::DRef &B, + const nb::DRef &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)"); +} diff --git a/tests/test_all.py b/tests/test_all.py index a3f5ba7d..b11bcb04 100644 --- a/tests/test_all.py +++ b/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)