diff --git a/src/connect_boundary_to_infinity.cpp b/src/connect_boundary_to_infinity.cpp new file mode 100644 index 00000000..f3b22f31 --- /dev/null +++ b/src/connect_boundary_to_infinity.cpp @@ -0,0 +1,57 @@ +#include "default_types.h" +#include +#include +#include +#include +#include + +namespace nb = nanobind; +using namespace nb::literals; + +namespace pyigl +{ + auto connect_boundary_to_infinity_F( + const nb::DRef &F) + { + Eigen::MatrixXI FO; + igl::connect_boundary_to_infinity(F, FO); + return FO; + } + + auto connect_boundary_to_infinity_VF( + const nb::DRef &V, + const nb::DRef &F) + { + Eigen::MatrixXN VO; + Eigen::MatrixXI FO; + igl::connect_boundary_to_infinity(V, F, VO, FO); + return std::make_tuple(VO, FO); + } +} + +void bind_connect_boundary_to_infinity(nb::module_ &m) +{ + m.def( + "connect_boundary_to_infinity", + &pyigl::connect_boundary_to_infinity_F, + "F"_a, +R"(Connect all boundary edges to a fictitious point at infinity, closing an +open mesh (the result is edge-manifold if the input was manifold). + +@param[in] F #F by 3 list of face indices into some V +@param[out] FO #F+#O by 3 list of face indices into [V; inf], original F come + first)"); + m.def( + "connect_boundary_to_infinity", + &pyigl::connect_boundary_to_infinity_VF, + "V"_a, + "F"_a, +R"(Connect all boundary edges to a fictitious point at infinity, also +returning the augmented vertex list. + +@param[in] V #V by 3 list of vertex positions +@param[in] F #F by 3 list of face indices into rows of V +@param[out] VO #V+1 by 3 list of vertex positions; original V come first and + the last row is (inf, inf, inf) +@param[out] FO #F+#O by 3 list of face indices into rows of VO)"); +} diff --git a/src/fit_rotations.cpp b/src/fit_rotations.cpp new file mode 100644 index 00000000..1006fe94 --- /dev/null +++ b/src/fit_rotations.cpp @@ -0,0 +1,35 @@ +#include "default_types.h" +#include +#include +#include +#include + +namespace nb = nanobind; +using namespace nb::literals; + +namespace pyigl +{ + auto fit_rotations( + const nb::DRef &S, + const bool single_precision) + { + Eigen::MatrixXN R; + igl::fit_rotations(S, single_precision, R); + return R; + } +} + +void bind_fit_rotations(nb::module_ &m) +{ + m.def( + "fit_rotations", + &pyigl::fit_rotations, + "S"_a, + "single_precision"_a = false, +R"(Given a stack of covariance matrices, find the closest rotation to each +(e.g. the local step of an as-rigid-as-possible optimization). + +@param[in] S (n*dim) by dim stack of covariance matrices +@param[in] single_precision whether to use single precision (faster) +@param[out] R dim by (dim*n) list of rotations)"); +} diff --git a/src/path_to_edges.cpp b/src/path_to_edges.cpp new file mode 100644 index 00000000..ff9f0544 --- /dev/null +++ b/src/path_to_edges.cpp @@ -0,0 +1,35 @@ +#include "default_types.h" +#include +#include +#include +#include + +namespace nb = nanobind; +using namespace nb::literals; + +namespace pyigl +{ + auto path_to_edges( + const nb::DRef &I, + const bool make_loop) + { + Eigen::MatrixXI E; + igl::path_to_edges(I, E, make_loop); + return E; + } +} + +void bind_path_to_edges(nb::module_ &m) +{ + m.def( + "path_to_edges", + &pyigl::path_to_edges, + "I"_a, + "make_loop"_a = false, +R"(Given an ordered list of vertex indices describing a path, return the list +of undirected edges connecting consecutive entries. + +@param[in] I #I list of vertex indices along the path +@param[in] make_loop whether to also connect the last vertex back to the first +@param[out] E #E by 2 list of edges (E = #I-1, or #I if make_loop))"); +} diff --git a/src/planarize_quad_mesh.cpp b/src/planarize_quad_mesh.cpp new file mode 100644 index 00000000..4cdf42ff --- /dev/null +++ b/src/planarize_quad_mesh.cpp @@ -0,0 +1,43 @@ +#include "default_types.h" +#include +#include +#include +#include + +namespace nb = nanobind; +using namespace nb::literals; + +namespace pyigl +{ + auto planarize_quad_mesh( + const nb::DRef &Vin, + const nb::DRef &F, + const int maxIter, + const Numeric threshold) + { + // Vout must share Vin's (deduced) type, so pass a concrete copy of the input + Eigen::MatrixXN Vin_copy = Vin; + Eigen::MatrixXN Vout; + igl::planarize_quad_mesh(Vin_copy, F, maxIter, threshold, Vout); + return Vout; + } +} + +void bind_planarize_quad_mesh(nb::module_ &m) +{ + m.def( + "planarize_quad_mesh", + &pyigl::planarize_quad_mesh, + "Vin"_a, + "F"_a, + "maxIter"_a = 100, + "threshold"_a = 0.005, +R"(Planarize a quad mesh, moving vertices so that each quad becomes as planar +as possible while staying close to the original surface. + +@param[in] Vin #V by 3 list of input mesh vertex positions +@param[in] F #F by 4 list of quad face indices into rows of Vin +@param[in] maxIter maximum number of optimization iterations +@param[in] threshold planarity threshold at which to stop +@param[out] Vout #V by 3 list of planarized vertex positions)"); +} diff --git a/src/point_simplex_squared_distance.cpp b/src/point_simplex_squared_distance.cpp new file mode 100644 index 00000000..dd0194e3 --- /dev/null +++ b/src/point_simplex_squared_distance.cpp @@ -0,0 +1,64 @@ +#include "default_types.h" +#include +#include +#include +#include +#include + +namespace nb = nanobind; +using namespace nb::literals; + +namespace pyigl +{ + auto point_simplex_squared_distance( + const nb::DRef &p, + const nb::DRef &V, + const nb::DRef &Ele, + const Integer i) + { + // Use concrete Eigen types (template deduction is finicky with DRef maps here). + // DIM (the ambient dimension) is a non-deduced template parameter, so we must + // supply it explicitly; dispatch on the number of columns of V. + const Eigen::RowVectorXN pc = p; + const Eigen::MatrixXN Vc = V; + const Eigen::MatrixXI Elec = Ele; + Numeric sqr_d; + Eigen::RowVectorXN c; + Eigen::RowVectorXN b; + switch (Vc.cols()) + { + case 2: + igl::point_simplex_squared_distance<2>(pc, Vc, Elec, i, sqr_d, c, b); + break; + case 3: + igl::point_simplex_squared_distance<3>(pc, Vc, Elec, i, sqr_d, c, b); + break; + default: + throw std::runtime_error( + "point_simplex_squared_distance: only 2D and 3D points are supported"); + } + return std::make_tuple(sqr_d, c, b); + } +} + +void bind_point_simplex_squared_distance(nb::module_ &m) +{ + m.def( + "point_simplex_squared_distance", + &pyigl::point_simplex_squared_distance, + "p"_a, + "V"_a, + "Ele"_a, + "i"_a, +R"(Determine the squared distance from a point p to the i-th linear simplex +(point/segment/triangle) of a mesh, along with the closest point and its +barycentric coordinates. + +@param[in] p dim-long query point +@param[in] V #V by dim list of simplex-corner vertex positions +@param[in] Ele #Ele by (1|2|3) list of simplex indices into rows of V +@param[in] i index of the simplex (row of Ele) to measure against +@param[out] sqr_d squared distance from p to the i-th simplex +@param[out] c dim-long closest point on the i-th simplex +@param[out] b barycentric coordinates of the closest point c)"); +} diff --git a/src/polar_svd3x3.cpp b/src/polar_svd3x3.cpp new file mode 100644 index 00000000..a77a37ed --- /dev/null +++ b/src/polar_svd3x3.cpp @@ -0,0 +1,32 @@ +#include "default_types.h" +#include +#include +#include +#include + +namespace nb = nanobind; +using namespace nb::literals; + +namespace pyigl +{ + auto polar_svd3x3(const nb::DRef &A) + { + Eigen::Matrix Am = A; + Eigen::Matrix R; + igl::polar_svd3x3(Am, R); + return Eigen::MatrixXN(R); + } +} + +void bind_polar_svd3x3(nb::module_ &m) +{ + m.def( + "polar_svd3x3", + &pyigl::polar_svd3x3, + "A"_a, +R"(Compute the closest rotation matrix R to a 3x3 matrix A via the polar +decomposition (using a fast fixed-size SVD). + +@param[in] A 3 by 3 matrix +@param[out] R 3 by 3 closest rotation matrix (det(R) = +1))"); +} diff --git a/src/quad_edges.cpp b/src/quad_edges.cpp new file mode 100644 index 00000000..0d74f7b6 --- /dev/null +++ b/src/quad_edges.cpp @@ -0,0 +1,33 @@ +#include "default_types.h" +#include +#include +#include +#include + +namespace nb = nanobind; +using namespace nb::literals; + +namespace pyigl +{ + auto quad_edges(const nb::DRef &Q) + { + // igl::quad_edges is instantiated only for 32-bit int matrices (it builds an + // Eigen::MatrixXi internally), so cast to int and back. + Eigen::Matrix Qi = Q.cast(); + Eigen::Matrix Ei; + igl::quad_edges(Qi, Ei); + return Eigen::MatrixXI(Ei.cast()); + } +} + +void bind_quad_edges(nb::module_ &m) +{ + m.def( + "quad_edges", + &pyigl::quad_edges, + "Q"_a, +R"(Compute the list of unique undirected edges of a quad mesh. + +@param[in] Q #Q by 4 list of quad face indices +@param[out] E #E by 2 list of unique edges)"); +} diff --git a/src/quad_planarity.cpp b/src/quad_planarity.cpp new file mode 100644 index 00000000..d4e8c157 --- /dev/null +++ b/src/quad_planarity.cpp @@ -0,0 +1,36 @@ +#include "default_types.h" +#include +#include +#include +#include + +namespace nb = nanobind; +using namespace nb::literals; + +namespace pyigl +{ + auto quad_planarity( + const nb::DRef &V, + const nb::DRef &F) + { + Eigen::VectorXN P; + igl::quad_planarity(V, F, P); + return P; + } +} + +void bind_quad_planarity(nb::module_ &m) +{ + m.def( + "quad_planarity", + &pyigl::quad_planarity, + "V"_a, + "F"_a, +R"(Compute a planarity value for each quad in a quad mesh. Planarity is the +distance between the two diagonals of the quad divided by the average of the +diagonals' lengths (0 for a perfectly planar quad). + +@param[in] V #V by 3 list of mesh vertex positions +@param[in] F #F by 4 list of quad face indices into rows of V +@param[out] P #F list of planarity values)"); +} diff --git a/src/ramer_douglas_peucker.cpp b/src/ramer_douglas_peucker.cpp new file mode 100644 index 00000000..637030ea --- /dev/null +++ b/src/ramer_douglas_peucker.cpp @@ -0,0 +1,38 @@ +#include "default_types.h" +#include +#include +#include +#include +#include + +namespace nb = nanobind; +using namespace nb::literals; + +namespace pyigl +{ + auto ramer_douglas_peucker( + const nb::DRef &P, + const Numeric tol) + { + Eigen::MatrixXN S; + Eigen::VectorXI J; + igl::ramer_douglas_peucker(P, tol, S, J); + return std::make_tuple(S, J); + } +} + +void bind_ramer_douglas_peucker(nb::module_ &m) +{ + m.def( + "ramer_douglas_peucker", + &pyigl::ramer_douglas_peucker, + "P"_a, + "tol"_a, +R"(Ramer-Douglas-Peucker piecewise-linear curve simplification. + +@param[in] P #P by dim ordered list of vertices along the curve +@param[in] tol tolerance (maximal Euclidean distance allowed between the new + line and a removed vertex) +@param[out] S #S by dim ordered list of points along the simplified curve +@param[out] J #S list of indices into P so that S = P[J,:])"); +} diff --git a/src/smooth_corner_adjacency.cpp b/src/smooth_corner_adjacency.cpp new file mode 100644 index 00000000..d391ab92 --- /dev/null +++ b/src/smooth_corner_adjacency.cpp @@ -0,0 +1,44 @@ +#include "default_types.h" +#include +#include +#include +#include +#include + +namespace nb = nanobind; +using namespace nb::literals; + +namespace pyigl +{ + auto smooth_corner_adjacency( + const nb::DRef &V, + const nb::DRef &F, + const Numeric corner_threshold_radians) + { + Eigen::VectorXI CI; + Eigen::VectorXI CC; + igl::smooth_corner_adjacency(V, F, corner_threshold_radians, CI, CC); + return std::make_tuple(CI, CC); + } +} + +void bind_smooth_corner_adjacency(nb::module_ &m) +{ + m.def( + "smooth_corner_adjacency", + &pyigl::smooth_corner_adjacency, + "V"_a, + "F"_a, + "corner_threshold_radians"_a, +R"(Determine the corner-to-face adjacency relationship that can be used for +computing crease-aware per-corner normals: two faces incident on a shared +vertex are grouped together for that corner only if their dihedral angle is +below the given threshold. + +@param[in] V #V by 3 list of mesh vertex positions +@param[in] F #F by 3 list of triangle indices into rows of V +@param[in] corner_threshold_radians dihedral angle threshold in radians +@param[out] CI flat list of face neighbors, indexed via CC (CSR-style) +@param[out] CC #F*3+1 cumulative sum of corner-neighbor counts so that the + neighbors of corner c are CI[CC[c]:CC[c+1]])"); +} diff --git a/tests/test_all.py b/tests/test_all.py index b11bcb04..ed893aa5 100644 --- a/tests/test_all.py +++ b/tests/test_all.py @@ -1429,6 +1429,133 @@ def test_new_triangle_algorithms(): assert V_ref.shape[0] >= V_tri.shape[0] +def test_polar_svd3x3_and_fit_rotations(): + # A = R * S with R a known rotation and S SPD; polar decomposition recovers R + theta = 0.7 + c, s = np.cos(theta), np.sin(theta) + R = np.array([[c, -s, 0.0], [s, c, 0.0], [0.0, 0.0, 1.0]]) + # Diagonal (SPD) stretch => the closest rotation to A = R @ S is exactly R + S = np.diag([2.0, 1.5, 1.0]) + A = R @ S + Rr = igl.polar_svd3x3(A) + assert Rr.shape == (3, 3) + # A proper rotation that recovers R + np.testing.assert_allclose(Rr @ Rr.T, np.eye(3), atol=1e-5) + assert np.isclose(np.linalg.det(Rr), 1.0) + np.testing.assert_allclose(Rr, R, atol=1e-5) + # Left polar factor A @ Rr.T is symmetric positive-(semi)definite + P = A @ Rr.T + np.testing.assert_allclose(P, P.T, atol=1e-5) + + # fit_rotations wraps polar_svd3x3 on a single covariance (dim by dim*n, n=1). + # It stores each rotation transposed relative to polar_svd3x3 (ARAP convention). + Rf = igl.fit_rotations(A) + assert Rf.shape == (3, 3) + np.testing.assert_allclose(Rf @ Rf.T, np.eye(3), atol=1e-5) + assert np.isclose(np.linalg.det(Rf), 1.0) + np.testing.assert_allclose(Rf, Rr.T, atol=1e-5) + + +def test_point_simplex_squared_distance(): + V = np.array([[0.0, 0.0, 0.0], + [1.0, 0.0, 0.0], + [0.0, 1.0, 0.0]], dtype=np.float64) + Ele = np.array([[0, 1, 2]], dtype=np.int64) + # Point 2 above the interior point (0.25,0.25,0) + p = np.array([0.25, 0.25, 2.0]) + sqr_d, c, b = igl.point_simplex_squared_distance(p, V, Ele, 0) + assert np.isclose(sqr_d, 4.0) + np.testing.assert_allclose(c, [0.25, 0.25, 0.0], atol=1e-9) + assert np.isclose(b.sum(), 1.0) + np.testing.assert_allclose(b @ V, c, atol=1e-9) + + # 2D case: point-to-segment + V2 = np.array([[0.0, 0.0], [2.0, 0.0]], dtype=np.float64) + E2 = np.array([[0, 1]], dtype=np.int64) + sqr2, c2, b2 = igl.point_simplex_squared_distance(np.array([1.0, 3.0]), V2, E2, 0) + assert np.isclose(sqr2, 9.0) + np.testing.assert_allclose(c2, [1.0, 0.0], atol=1e-9) + + +def test_quad_mesh_helpers(): + # Unit planar square as a single quad + V = np.array([[0.0, 0.0, 0.0], + [1.0, 0.0, 0.0], + [1.0, 1.0, 0.0], + [0.0, 1.0, 0.0]], dtype=np.float64) + F = np.array([[0, 1, 2, 3]], dtype=np.int64) + + # quad_edges: 4 unique edges of the quad + E = igl.quad_edges(F) + assert E.shape == (4, 2) + edge_set = set(tuple(sorted(e)) for e in E.tolist()) + assert edge_set == {(0, 1), (1, 2), (2, 3), (0, 3)} + + # quad_planarity: a planar quad has ~zero non-planarity + P = igl.quad_planarity(V, F) + assert P.shape[0] == 1 + assert np.isclose(P[0], 0.0, atol=1e-9) + + # planarize_quad_mesh: a non-planar quad becomes more planar + Vnp = V.copy() + Vnp[2, 2] = 0.5 # lift one corner out of plane + P_before = igl.quad_planarity(Vnp, F)[0] + Vout = igl.planarize_quad_mesh(Vnp, F, maxIter=100, threshold=1e-4) + assert Vout.shape == V.shape + P_after = igl.quad_planarity(Vout, F)[0] + assert P_after <= P_before + 1e-12 + + +def test_ramer_douglas_peucker(): + # Nearly-collinear first three points collapse; the sharp turn is kept + P = np.array([[0.0, 0.0], + [1.0, 0.0001], + [2.0, 0.0], + [2.0, 2.0]], dtype=np.float64) + S, J = igl.ramer_douglas_peucker(P, 0.01) + # Interior near-collinear vertex (index 1) is dropped + assert J.tolist() == [0, 2, 3] + np.testing.assert_allclose(S, P[J], atol=1e-12) + + +def test_path_to_edges_and_connect_boundary(): + # path_to_edges: open path and closed loop + I = np.array([0, 1, 2, 3], dtype=np.int64) + E = igl.path_to_edges(I) + assert E.tolist() == [[0, 1], [1, 2], [2, 3]] + E_loop = igl.path_to_edges(np.array([0, 1, 2], dtype=np.int64), make_loop=True) + assert E_loop.tolist() == [[0, 1], [1, 2], [2, 0]] + + # connect_boundary_to_infinity on an open mesh (square from two triangles) + V = np.array([[0.0, 0.0, 0.0], + [1.0, 0.0, 0.0], + [1.0, 1.0, 0.0], + [0.0, 1.0, 0.0]], dtype=np.float64) + F = np.array([[0, 1, 2], [0, 2, 3]], dtype=np.int64) + # (V,F) overload appends a single new vertex at infinity + VO, FO = igl.connect_boundary_to_infinity(V, F) + assert VO.shape[0] == V.shape[0] + 1 + assert FO.shape[0] == F.shape[0] + 4 # 4 boundary edges connected + np.testing.assert_allclose(VO[:V.shape[0]], V) + # F-only overload returns the same number of faces + FO2 = igl.connect_boundary_to_infinity(F) + assert FO2.shape[0] == F.shape[0] + 4 + + +def test_smooth_corner_adjacency(): + # Two triangles sharing an edge; CSR-style corner adjacency + V = np.array([[0.0, 0.0, 0.0], + [1.0, 0.0, 0.0], + [1.0, 1.0, 0.0], + [0.0, 1.0, 0.0]], dtype=np.float64) + F = np.array([[0, 1, 2], [0, 2, 3]], dtype=np.int64) + CI, CC = igl.smooth_corner_adjacency(V, F, 0.35) + # CC is a cumulative-sum offset array of length 3*#F + 1 + assert CC.shape[0] == 3 * F.shape[0] + 1 + assert CC[0] == 0 + assert CC[-1] == CI.shape[0] + assert np.all(np.diff(CC) >= 0) + def test_remesh_at_points(): # Single triangle in 3D V = np.array([[0.0, 0.0, 0.0],