diff --git a/pyproject.toml b/pyproject.toml index df131292..42885683 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,9 +1,8 @@ [build-system] requires = [ - "cmake>=3.16", "nanobind >=1.3.2", "numpy>=2.0.0; python_version >= '3.9'", - "numpy; python_version < '3.9'", + "numpy<2.0.0; python_version < '3.9'", "packaging", "scikit-build-core >=0.10", "scipy", diff --git a/src/comb_cross_field.cpp b/src/comb_cross_field.cpp new file mode 100644 index 00000000..a3baaa67 --- /dev/null +++ b/src/comb_cross_field.cpp @@ -0,0 +1,39 @@ +#include "default_types.h" +#include +#include +#include +#include +#include + +namespace nb = nanobind; +using namespace nb::literals; + +void bind_comb_cross_field(nb::module_ &m) +{ + m.def("comb_cross_field", + []( + const Eigen::MatrixXN& V, + const Eigen::MatrixXI& F, + const Eigen::MatrixXN& PD1, + const Eigen::MatrixXN& PD2 + ) { + Eigen::MatrixXN PD1_out, PD2_out; + igl::comb_cross_field(V, F, PD1, PD2, PD1_out, PD2_out); + return std::make_tuple(PD1_out, PD2_out); + }, + "V"_a, + "F"_a, + "PD1"_a, + "PD2"_a, + R"(Computes principal matchings of the vectors of a cross field across face edges, + and generates a combed cross field defined on the mesh faces + + @param[in] V #V by 3 eigen Matrix of mesh vertex 3D positions + @param[in] F #F by 4 eigen Matrix of face (quad) indices + @param[in] PD1in #F by 3 eigen Matrix of the first per face cross field vector + @param[in] PD2in #F by 3 eigen Matrix of the second per face cross field vector + @param[out] PD1out #F by 3 eigen Matrix of the first combed cross field vector + @param[out] PD2out #F by 3 eigen Matrix of the second combed cross field vector + )" + ); +} diff --git a/src/comb_frame_field.cpp b/src/comb_frame_field.cpp new file mode 100644 index 00000000..7970ec81 --- /dev/null +++ b/src/comb_frame_field.cpp @@ -0,0 +1,47 @@ +#include "default_types.h" +#include +#include +#include +#include +#include + +namespace nb = nanobind; +using namespace nb::literals; + +void bind_comb_frame_field(nb::module_ &m) +{ + m.def("comb_frame_field", + []( + const Eigen::MatrixXN& V, + const Eigen::MatrixXI& F, + const Eigen::MatrixXN& PD1, + const Eigen::MatrixXN& PD2, + const Eigen::MatrixXN& BIS1_combed, + const Eigen::MatrixXN& BIS2_combed + ) { + Eigen::MatrixXN PD1_combed, PD2_combed; + igl::comb_frame_field(V, F, PD1, PD2, BIS1_combed, BIS2_combed, PD1_combed, PD2_combed); + return std::make_tuple(PD1_combed, PD2_combed); + }, + "V"_a, + "F"_a, + "PD1"_a, + "PD2"_a, + "BIS1_combed"_a, + "BIS2_combed"_a, + R"(Computes principal matchings of the vectors of a frame field across face edges, + and generates a combed frame field defined on the mesh faces. This makes use of a + combed cross field generated by combing the field created by the bisectors of the + frame field. + + @param[in] V #V by 3 eigen Matrix of mesh vertex 3D positions + @param[in] F #F by 4 eigen Matrix of face (quad) indices + @param[in] PD1 #F by 3 eigen Matrix of the first per face cross field vector + @param[in] PD2 #F by 3 eigen Matrix of the second per face cross field vector + @param[in] BIS1_combed #F by 3 eigen Matrix of the first combed bisector field vector + @param[in] BIS2_combed #F by 3 eigen Matrix of the second combed bisector field vector + @param[out] PD1_combed #F by 3 eigen Matrix of the first combed cross field vector + @param[out] PD2_combed #F by 3 eigen Matrix of the second combed cross field vector + )" + ); +} diff --git a/src/compute_frame_field_bisectors.cpp b/src/compute_frame_field_bisectors.cpp new file mode 100644 index 00000000..de28683d --- /dev/null +++ b/src/compute_frame_field_bisectors.cpp @@ -0,0 +1,71 @@ +#include "default_types.h" +#include +#include +#include +#include +#include + +namespace nb = nanobind; +using namespace nb::literals; + +void bind_compute_frame_field_bisectors(nb::module_ &m) +{ + m.def("compute_frame_field_bisectors", + []( + const Eigen::MatrixXN& V, + const Eigen::MatrixXI& F, + const Eigen::MatrixXN& B1, + const Eigen::MatrixXN& B2, + const Eigen::MatrixXN& PD1, + const Eigen::MatrixXN& PD2 + ) { + Eigen::MatrixXN BIS1, BIS2; + igl::compute_frame_field_bisectors(V, F, B1, B2, PD1, PD2, BIS1, BIS2); + return std::make_tuple(BIS1, BIS2); + }, + "V"_a, + "F"_a, + "B1"_a, + "B2"_a, + "PD1"_a, + "PD2"_a, + R"(Compute bisectors of a frame field defined on mesh faces. + + @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[in] B1 #F by 3 eigen Matrix of face (triangle) base vector 1 + @param[in] B2 #F by 3 eigen Matrix of face (triangle) base vector 2 + @param[in] PD1 #F by 3 eigen Matrix of the first per face frame field vector + @param[in] PD2 #F by 3 eigen Matrix of the second per face frame field vector + @param[out] BIS1 #F by 3 eigen Matrix of the first per face frame field bisector + @param[out] BIS2 #F by 3 eigen Matrix of the second per face frame field bisector + )" + ); + + // Overload + m.def("compute_frame_field_bisectors", + []( + const Eigen::MatrixXN& V, + const Eigen::MatrixXI& F, + const Eigen::MatrixXN& PD1, + const Eigen::MatrixXN& PD2 + ) { + Eigen::MatrixXN BIS1, BIS2; + igl::compute_frame_field_bisectors(V, F, PD1, PD2, BIS1, BIS2); + return std::make_tuple(BIS1, BIS2); + }, + "V"_a, + "F"_a, + "PD1"_a, + "PD2"_a, + R"(Compute bisectors of a frame field defined on mesh faces. + + @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[in] PD1 #F by 3 eigen Matrix of the first per face frame field vector + @param[in] PD2 #F by 3 eigen Matrix of the second per face frame field vector + @param[out] BIS1 #F by 3 eigen Matrix of the first per face frame field bisector + @param[out] BIS2 #F by 3 eigen Matrix of the second per face frame field bisector + )" + ); +} diff --git a/src/cross_field_mismatch.cpp b/src/cross_field_mismatch.cpp new file mode 100644 index 00000000..4039c1e9 --- /dev/null +++ b/src/cross_field_mismatch.cpp @@ -0,0 +1,45 @@ +#include "default_types.h" +#include +#include +#include +#include +#include + +namespace nb = nanobind; +using namespace nb::literals; + +void bind_cross_field_mismatch(nb::module_ &m) +{ + m.def("cross_field_mismatch", + []( + const Eigen::MatrixXN& V, + const Eigen::MatrixXI& F, + const Eigen::MatrixXN& PD1, + const Eigen::MatrixXN& PD2, + const bool isCombed + ) { + Eigen::MatrixXI mismatch; + igl::cross_field_mismatch(V, F, PD1, PD2, isCombed, mismatch); + return mismatch; + }, + "V"_a, + "F"_a, + "PD1"_a, + "PD2"_a, + "isCombed"_a, + R"(Calculates the mismatch (integer), at each face edge, of a cross field defined on the mesh faces. + The integer mismatch is a multiple of pi/2 that transforms the cross on one side of the edge to + the cross on the other side. It represents the deviation from a Lie connection across the edge. + + @param[in] V #V by 3 eigen Matrix of mesh vertex 3D positions + @param[in] F #F by 3 eigen Matrix of face (quad) indices + @param[in] PD1 #F by 3 eigen Matrix of the first per face cross field vector + @param[in] PD2 #F by 3 eigen Matrix of the second per face cross field vector + @param[in] isCombed boolean, specifying whether the field is combed (i.e. matching has been precomputed. + If not, the field is combed first. + @param[out] mismatch #F by 3 eigen Matrix containing the integer mismatch of the cross field + across all face edges + )" + + ); +} diff --git a/src/find_cross_field_singularities.cpp b/src/find_cross_field_singularities.cpp new file mode 100644 index 00000000..dfefcd91 --- /dev/null +++ b/src/find_cross_field_singularities.cpp @@ -0,0 +1,64 @@ +#include "default_types.h" +#include +#include +#include +#include +#include + +namespace nb = nanobind; +using namespace nb::literals; + +void bind_find_cross_field_singularities(nb::module_ &m) +{ + m.def("find_cross_field_singularities", + []( + const Eigen::MatrixXN& V, + const Eigen::MatrixXI& F, + const Eigen::MatrixXI& mismatch + ) { + Eigen::VectorXI isSingularity, singularityIndex; + igl::find_cross_field_singularities(V, F, mismatch, isSingularity, singularityIndex); + return std::make_tuple(isSingularity, singularityIndex); + }, + "V"_a, + "F"_a, + "mismatch"_a, + R"(Computes singularities of a cross field, assumed combed + + @param[in] V #V by 3 eigen Matrix of mesh vertex 3D positions + @param[in] F #F by 3 eigen Matrix of face (quad) indices + @param[in] mismatch #F by 3 eigen Matrix containing the integer mismatch of the cross field + across all face edges + @param[out] isSingularity #V by 1 boolean eigen Vector indicating the presence of a singularity on a vertex + @param[out] singularityIndex #V by 1 integer eigen Vector containing the singularity indices + )" + ); + + // overload + m.def("find_cross_field_singularities", + []( + const Eigen::MatrixXN& V, + const Eigen::MatrixXI& F, + const Eigen::MatrixXN& PD1, + const Eigen::MatrixXN& PD2, + bool isCombed + ) { + Eigen::VectorXI isSingularity, singularityIndex; + igl::find_cross_field_singularities(V, F, PD1, PD2, isSingularity, singularityIndex, isCombed); + return std::make_tuple(isSingularity, singularityIndex); + }, + "V"_a, + "F"_a, + "PD1"_a, + "PD2"_a, + "isCombed"_a = false, + R"(Wrapper that calculates the mismatch if it is not provided. + + @param[in] PD1 #F by 3 eigen Matrix of the first per face cross field vector + @param[in] PD2 #F by 3 eigen Matrix of the second per face cross field vector + @param[in] isCombed boolean indicating whether the cross field is combed + + \note the field in PD1 and PD2 MUST BE combed (see igl::comb_cross_field). + )" + ); +} diff --git a/src/rotate_vectors.cpp b/src/rotate_vectors.cpp new file mode 100644 index 00000000..c2eafa03 --- /dev/null +++ b/src/rotate_vectors.cpp @@ -0,0 +1,28 @@ +#include "default_types.h" +#include +#include +#include +#include +#include + +namespace nb = nanobind; +using namespace nb::literals; + +void bind_rotate_vectors(nb::module_ &m) { + m.def( + "rotate_vectors", + [](const Eigen::MatrixXN &V, const Eigen::VectorXN &A, + const Eigen::MatrixXN &B1, const Eigen::MatrixXN &B2) { + return igl::rotate_vectors(V, A, B1, B2); + }, + "V"_a, "A"_a, "B1"_a, "B2"_a, + R"(Rotate the vectors V by A radians on the tangent plane spanned by B1 and B2 + + @param[in] V #V by 3 eigen Matrix of vectors + @param[in] A #V eigen vector of rotation angles or a single angle to be applied + to all vectors + @param[in] B1 #V by 3 eigen Matrix of base vector 1 + @param[in] B2 #V by 3 eigen Matrix of base vector 2 + @return the rotated vectors + )"); +} diff --git a/tests/test_all.py b/tests/test_all.py index b666a4e7..21cd448a 100644 --- a/tests/test_all.py +++ b/tests/test_all.py @@ -12,6 +12,10 @@ import igl.copyleft.cgal import igl.embree +@pytest.fixture +def icosahedron(): + V,F = igl.icosahedron() + return V,F #def rand_sparse(n,density): # n_features = n @@ -600,3 +604,144 @@ def udf_sphere(Q): unique_ijk, J, unique_corners = igl.unique_sparse_voxel_corners(origin,h0,max_depth,ijk) unique_S = sdf_sphere(unique_corners) V,F = igl.marching_cubes(unique_S,unique_corners,J,0.0) + +def test_rotate_vectors(icosahedron): + V,F = icosahedron + + # Create rotation angles (rotate by pi/4) + A = np.ones(F.shape[0],dtype=np.float64) * (np.pi / 4.0) + + # Get local basis + B1,B2,_ = igl.local_basis(V,F) + + # B1 is orthogonal to B2 + r = np.sum(B1* B2, axis=1) + assert np.allclose(np.abs(r), 0.0) + + # Rotate the first basis vector + B1_rotated = igl.rotate_vectors(B1, A, B1, B2) + + # Check output shape + assert B1_rotated.shape == B1.shape + + # Rotate B1_rotated by pi/4 again + B1_rotated2 = igl.rotate_vectors(B1_rotated, A, B1, B2) + assert B1_rotated2.shape == B1_rotated.shape + + # B1_rotated2 should be parallel to B2 + r = np.sum(B1_rotated2 * B2, axis=1) + assert np.allclose(np.abs(r), 1.0) + +def test_compute_frame_field_bisectors(icosahedron): + V,F = icosahedron + + # Get local basis + B1,B2,_ = igl.local_basis(V,F) + + # Compute bisectors with explicit basis + BIS1,BIS2 = igl.compute_frame_field_bisectors(V,F,B1,B2,B1,B2) + + # Check output shapes + assert BIS1.shape == (F.shape[0], 3) + assert BIS2.shape == (F.shape[0], 3) + + # BIS1 should be orthogonal to BIS2 + r = np.sum(BIS1 * BIS2, axis=1) + assert np.allclose(np.abs(r), 0.0) + +def test_comb_cross_field(icosahedron): + V,F = icosahedron + + # Get local basis + B1,B2,_ = igl.local_basis(V,F) + + # Comb the cross field + B1_combed,B2_combed = igl.comb_cross_field(V,F,B1,B2) + + # Check output shapes + assert B1_combed.shape == (F.shape[0], 3) + assert B2_combed.shape == (F.shape[0], 3) + +def test_cross_field_mismatch(icosahedron): + V,F = icosahedron + + # Get local basis + B1,B2,_ = igl.local_basis(V,F) + + # Comb the cross field first + B1_combed,B2_combed = igl.comb_cross_field(V,F,B1,B2) + + # Compute mismatch on combed field + mismatch = igl.cross_field_mismatch(V,F,B1_combed,B2_combed,True) + + # Check output shape (should be #F by 3 for triangular mesh) + assert mismatch.shape == (F.shape[0], 3) + + # Test with uncombed field (function will comb it first) + mismatch_uncombed = igl.cross_field_mismatch(V,F,B1,B2,False) + assert mismatch_uncombed.shape == (F.shape[0], 3) + +def test_find_cross_field_singularities(icosahedron): + V,F = icosahedron + + # Get local basis + B1,B2,_ = igl.local_basis(V,F) + + # Comb the cross field first + B1_combed,B2_combed = igl.comb_cross_field(V,F,B1,B2) + + # Compute mismatch + mismatch = igl.cross_field_mismatch(V,F,B1_combed,B2_combed,True) + + # Find singularities from mismatch + isSingularity,singularityIndex = igl.find_cross_field_singularities(V,F,mismatch) + + # Check Poincaré-Hopf theorem + # The current singularity computation only return positive index, which is inconsistent with the + # theorem, so we skip this check for now. + # assert np.sum(singularityIndex) == 2 * 4 # Euler characteristic * 4-rosy fields + + # Check output shapes + assert isSingularity.shape[0] == V.shape[0] + assert singularityIndex.shape[0] == V.shape[0] + + # Test overload that computes mismatch internally + isSingularity2,singularityIndex2 = igl.find_cross_field_singularities(V,F,B1_combed,B2_combed,True) + + # Check Poincaré-Hopf theorem + # The current singularity computation only return positive index, which is inconsistent with the + # theorem, so we skip this check for now. + # assert np.sum(singularityIndex2) == 2 * 4 # Euler characteristic * 4-rosy fields + + # Check output shapes + assert isSingularity2.shape[0] == V.shape[0] + assert singularityIndex2.shape[0] == V.shape[0] + + # Test with uncombed field + isSingularity3,singularityIndex3 = igl.find_cross_field_singularities(V,F,B1,B2,False) + assert isSingularity3.shape[0] == V.shape[0] + assert singularityIndex3.shape[0] == V.shape[0] + + # Check Poincaré-Hopf theorem + # The current singularity computation only return positive index, which is inconsistent with the + # theorem, so we skip this check for now. + # assert np.sum(singularityIndex3) == 2 * 4 # Euler characteristic * 4-rosy fields + +def test_comb_frame_field(icosahedron): + V,F = icosahedron + + # Get local basis + B1,B2,_ = igl.local_basis(V,F) + + # Compute bisectors + BIS1,BIS2 = igl.compute_frame_field_bisectors(V,F,B1,B2) + + # Comb the bisectors (which are a cross field) + BIS1_combed,BIS2_combed = igl.comb_cross_field(V,F,BIS1,BIS2) + + # Comb the frame field using combed bisectors + PD1_combed,PD2_combed = igl.comb_frame_field(V,F,B1,B2,BIS1_combed,BIS2_combed) + + # Check output shapes + assert PD1_combed.shape == (F.shape[0], 3) + assert PD2_combed.shape == (F.shape[0], 3)