Merge pull request #239 from libigl/alecjacobson/bump-libigl-e7a6ed

bump libigl and update decimates

.github/workflows/wheels.yml

@@ -37,7 +37,7 @@ jobs:
       CIBW_TEST_COMMAND: "python {project}/tests/test_basic.py"
       CIBW_BUILD: "${{ matrix.cpversion }}-${{ matrix.os.cibw-arch }}"
       CIBW_TEST_SKIP: "*-macosx_arm64"
-      CIBW_ENVIRONMENT: "MAX_JOBS=${{ matrix.os.runs-on == 'macos-latest' && 3 || 2 }}"
+      CIBW_ENVIRONMENT: "MAX_JOBS=${{ matrix.os.runs-on == 'macos-latest' && 3 || 2 }} PIP_CONSTRAINT=constraints.txt"
       # Why universal2 here? It's not included above in CIBW_BUILD
       CIBW_ARCHS_MACOS: "x86_64 arm64 universal2"
       CIBW_ENVIRONMENT_MACOS: "MACOSX_DEPLOYMENT_TARGET=10.13 CMAKE_OSX_ARCHITECTURES=\"${{ matrix.os.cibw-arch == 'macosx_x86_64' && 'x86_64' || matrix.os.cibw-arch == 'macosx_arm64' && 'arm64' || matrix.os.cibw-arch == 'macosx_universal2' && 'arm64;x86_64' || '' }}\""
@@ -89,8 +89,9 @@ jobs:
           python -m pip install cibuildwheel
 
       - name: Build wheels
+        shell: bash
         run: |
-          python -m cibuildwheel --output-dir wheelhouse
+          PIP_CONSTRAINT=$GITHUB_WORKSPACE/constraints.txt python -m cibuildwheel --output-dir wheelhouse
 
       # Upload binaries to github
       - uses: actions/upload-artifact@v4

classes/AABB.cpp

@@ -4,58 +4,108 @@
 #include <pybind11/stl.h>
 #include <variant>
 #include <variant>
+#include <iostream>
 
 
 namespace py = pybind11;
+
+template <typename T> 
+void check_c_contiguity(py::array_t<T> V, const std::string name = "arg")
+{
+  auto V_buf = V.request();
+  if (V_buf.ndim != 2 || V_buf.strides[1] != sizeof(T)) {
+      throw std::runtime_error(name+" must be C-contiguous with a 2D shape.");
+  }
+};
+template <typename T> 
+void check_f_contiguity(py::array_t<T> V, const std::string name = "arg")
+{
+  auto V_buf = V.request();
+  if (V_buf.ndim != 2 || V_buf.strides[1] == sizeof(T)) {
+      throw std::runtime_error(name+" must be F-contiguous (ColMajor) with a 2D shape.");
+  }
+};
+
 template <int DIM>
 void init_AABB(py::module_ &m)
 {
-  using AABB_f64_DIM = igl::AABB<Eigen::MatrixXd,DIM>;
+  using namespace Eigen;
+  using MatrixXdRowMajor = Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor>;
+  using MatrixXiRowMajor = Eigen::Matrix<int,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor>;
+  using MapMatrixXdRowMajor = Eigen::Map<const MatrixXdRowMajor>;
+  using MapMatrixXiRowMajor = Eigen::Map<const MatrixXiRowMajor>;
+  using MapMatrixXd = Eigen::Map<const Eigen::MatrixXd>;
+  using MapMatrixXi = Eigen::Map<const Eigen::MatrixXi>;
+  using DerivedV = MapMatrixXd;
+  using AABB_f64_DIM = igl::AABB<DerivedV,DIM>;
   py::class_<AABB_f64_DIM >(m, (std::string("AABB_f64_")+std::to_string(DIM)).c_str() )
     .def(py::init([]()
       {
         std::unique_ptr<AABB_f64_DIM> self = std::make_unique<AABB_f64_DIM>();
         return self;
       }))
-    .def("init",[](AABB_f64_DIM & self, const Eigen::MatrixXd & V, const Eigen::MatrixXi & F)
+    .def("init", [](AABB_f64_DIM &self, py::array_t<double> V, py::array_t<int> F) 
       {
-        self.init(V,F);
-      },
-      py::arg("V"), py::arg("F"))
-    .def("squared_distance",[](
-          AABB_f64_DIM & self, 
-          const Eigen::MatrixXd & V, 
-          const Eigen::MatrixXi & F,
-          const Eigen::MatrixXd & P,
-          const bool return_index = false,
-          const bool return_closest_point = false) ->
-        std::variant<py::object,std::list<py::object> >
-      {
-        Eigen::VectorXd sqrD; 
-        Eigen::VectorXi I;
-        Eigen::MatrixXd C;
-        self.squared_distance(V,F,P,sqrD,I,C);
-        if(return_index && return_closest_point)
-        {
-          return std::list<py::object>({npe::move(sqrD),npe::move(I),npe::move(C)});
-        }else if(return_index)
-        {
-          return std::list<py::object>({npe::move(sqrD),npe::move(I)});
-        }else if(return_closest_point)
-        {
-          return std::list<py::object>({npe::move(sqrD),npe::move(C)});
-        }else
-        {
-          return npe::move(sqrD);
+        auto V_buf = V.request();
+        auto F_buf = F.request();
+
+        if (V_buf.ndim != 2 || F_buf.ndim != 2) {
+          throw std::runtime_error("Input matrices must be 2-dimensional.");
         }
-      },
-      py::arg("V"), 
-      py::arg("F"), 
-      py::arg("P"),
-      py::arg("return_index")=false,
-      py::arg("return_closest_point")=false
-      )
-    ;
+
+        check_f_contiguity(V,"V");
+        check_f_contiguity(F,"F");
+
+        MapMatrixXd V_eigen(static_cast<double*>(V_buf.ptr), V_buf.shape[0], V_buf.shape[1]);
+        MapMatrixXi F_eigen(static_cast<int*>(F_buf.ptr), F_buf.shape[0], F_buf.shape[1]);
+
+        self.init(V_eigen, F_eigen);
+      }, py::arg("V"), py::arg("F"))
+  .def("squared_distance",[](
+        AABB_f64_DIM & self, 
+        py::array_t<double> V, 
+        py::array_t<int> F,
+        py::array_t<double> P,
+        const bool return_index = false,
+        const bool return_closest_point = false) ->
+      std::variant<py::object,std::list<py::object> >
+    {
+      check_f_contiguity(V,"V");
+      check_f_contiguity(F,"F");
+      check_f_contiguity(P,"P");
+
+      auto V_buf = V.request();
+      auto F_buf = F.request();
+      auto P_buf = P.request();
+      MapMatrixXd V_eigen(static_cast<double*>(V_buf.ptr), V_buf.shape[0], V_buf.shape[1]);
+      MapMatrixXi F_eigen(static_cast<int*>(F_buf.ptr), F_buf.shape[0], F_buf.shape[1]);
+      MapMatrixXd P_eigen(static_cast<double*>(P_buf.ptr), P_buf.shape[0], P_buf.shape[1]);
+
+      Eigen::VectorXd sqrD; 
+      Eigen::VectorXi I;
+      MatrixXd C;
+      self.squared_distance(V_eigen,F_eigen,P_eigen,sqrD,I,C);
+      if(return_index && return_closest_point)
+      {
+        return std::list<py::object>({npe::move(sqrD),npe::move(I),npe::move(C)});
+      }else if(return_index)
+      {
+        return std::list<py::object>({npe::move(sqrD),npe::move(I)});
+      }else if(return_closest_point)
+      {
+        return std::list<py::object>({npe::move(sqrD),npe::move(C)});
+      }else
+      {
+        return npe::move(sqrD);
+      }
+    },
+    py::arg("V"), 
+    py::arg("F"), 
+    py::arg("P"),
+    py::arg("return_index")=false,
+    py::arg("return_closest_point")=false
+    )
+  ;
 }
 
 template void init_AABB<2>(py::module_ &);

cmake/PyiglDependencies.cmake

@@ -19,7 +19,7 @@ include(FetchContent)
 FetchContent_Declare(
     libigl
     GIT_REPOSITORY https://github.com/libigl/libigl.git
-    GIT_TAG v2.5.0
+    GIT_TAG f962e4a6b68afe978dc12a63702b7846a3e7a6ed
 )
 FetchContent_GetProperties(libigl)
 FetchContent_MakeAvailable(libigl)

constraints.txt

@@ -0,0 +1,2 @@
+numpy<=1.26.4
+

igl/_version.py

@@ -1,8 +1,8 @@
 # This file is part of libigl, a simple c++ geometry processing library.
 #
-# Copyright (C) 2023 Alec Jacobson
+# Copyright (C) 2024 Alec Jacobson
 #
 # This Source Code Form is subject to the terms of the Mozilla Public License
 # v. 2.0. If a copy of the MPL was not distributed with this file, You can
 # obtain one at http://mozilla.org/MPL/2.0/.
-__version__ = "2.5.1"
+__version__ = "2.5.4dev"

src/decimate.cpp

@@ -29,6 +29,7 @@ Parameters
 V  #V by dim list of vertex positions
 F  #F by 3 list of face indices into V.
 max_m  desired number of output faces
+block_intersections  whether to block intersections
 
 Returns
 -------
@@ -57,6 +58,7 @@ npe_doc(ds_decimate)
 npe_arg(v, dense_float, dense_double)
 npe_arg(f, dense_int32, dense_int64)
 npe_arg(max_m, size_t)
+npe_arg(block_intersections, bool)
 
 
 npe_begin_code()
@@ -69,7 +71,7 @@ npe_begin_code()
   Eigen::MatrixXi g;
   Eigen::VectorXi j;
   Eigen::VectorXi i;
-  bool reach = igl::decimate(v_copy, f_copy, max_m, u, g, j, i);
+  bool reach = igl::decimate(v_copy, f_copy, max_m, block_intersections, u, g, j, i);
   EigenDenseFloat u_row_major = u;
   EigenDenseInt g_row_major = g.template cast<typename EigenDenseInt::Scalar>();
   // FIXME: vector not allowing row major, but they should be essentially the same so i feel we can leave it as col major

src/in_element.cpp

@@ -11,7 +11,8 @@
 #include <igl/in_element.h>
 #include <igl/AABB.h>
 
-using AABB_f64_3 = igl::AABB<Eigen::MatrixXd,3>;
+using AABB_f64_3 = igl::AABB<Eigen::Map<const Eigen::MatrixXd>,3>;
+using AABB_f64_2 = igl::AABB<Eigen::Map<const Eigen::MatrixXd>,2>;
 const char* ds_in_element = R"igl_Qu8mg5v7(
 Determine whether each point in a list of points is in the elements of a mesh.
  
@@ -36,13 +37,13 @@ npe_arg(aabb, AABB_f64_3)
 npe_begin_code()
 
   Eigen::VectorXi I;
-  igl::in_element(V,Ele,Q,aabb,I);
+  Eigen::Map<const Eigen::MatrixXd> V_map(V.data(),V.rows(),V.cols());
+  igl::in_element(V_map,Ele,Q,aabb,I);
   return npe::move(I);
 
 npe_end_code()
 
 
-using AABB_f64_2 = igl::AABB<Eigen::MatrixXd,2>;
 npe_function(in_element_2)
 npe_doc(  ds_in_element)
 npe_arg(V, Eigen::MatrixXd)
@@ -52,7 +53,8 @@ npe_arg(aabb, AABB_f64_2)
 npe_begin_code()
 
   Eigen::VectorXi I;
-  igl::in_element(V,Ele,Q,aabb,I);
+  Eigen::Map<const Eigen::MatrixXd> V_map(V.data(),V.rows(),V.cols());
+  igl::in_element(V_map,Ele,Q,aabb,I);
   return npe::move(I);
 
 npe_end_code()

src/qslim.cpp

@@ -30,6 +30,7 @@ Parameters
 V  #V by dim list of vertex positions. Assumes that vertices w
 F  #F by 3 list of triangle indices into V
 max_m  desired number of output faces
+block_intersections  whether to block intersections
 
 
 Returns
@@ -60,6 +61,7 @@ npe_doc(ds_qslim)
 npe_arg(v, dense_float, dense_double)
 npe_arg(f, dense_int32, dense_int64)
 npe_arg(max_m, size_t)
+npe_arg(block_intersections, bool)
 
 
 npe_begin_code()
@@ -72,7 +74,7 @@ npe_begin_code()
   Eigen::MatrixXi g;
   Eigen::VectorXi j;
   Eigen::VectorXi i;
-  bool success = igl::qslim(v_copy, f_copy, max_m, u, g, j, i);
+  bool success = igl::qslim(v_copy, f_copy, max_m, block_intersections, u, g, j, i);
   EigenDenseFloat u_row_major = u;
   EigenDenseInt g_row_major = g.template cast<typename EigenDenseInt::Scalar>();
   // FIXME: vector not allowing row major, but they should be essentially the same so i feel we can leave it as col major

tests/test_basic.py

@@ -13,6 +13,8 @@
 import math
 import sys
 from git import Repo
+import faulthandler
+faulthandler.enable()
 
 
 
@@ -705,7 +707,7 @@ def test_cylinder(self):
         self.assertTrue(f.flags.c_contiguous)
 
     def test_decimate(self):
-        success, u, g, j, i = igl.decimate(self.v1, self.f1, 100)
+        success, u, g, j, i = igl.decimate(self.v1, self.f1, 100, False)
         self.assertEqual(u.shape[1], self.v1.shape[1])
         self.assertEqual(g.shape[1], 3)
         self.assertEqual(j.shape[0], g.shape[0])
@@ -823,7 +825,7 @@ def test_procrustes(self):
         self.assertTrue(r.dtype == t.dtype == self.v1.dtype)
 
     def test_qslim(self):
-        success, u, g, j, i = igl.qslim(self.v1, self.f1, 100)
+        success, u, g, j, i = igl.qslim(self.v1, self.f1, 100, False)
         self.assertEqual(u.dtype, self.v1.dtype)
         self.assertTrue(g.dtype == j.dtype == i.dtype == self.f1.dtype)
         self.assertEqual(u.shape[1], self.v1.shape[1])
@@ -2430,7 +2432,7 @@ def test_fast_winding_number_for_meshes(self):
 
     def test_flip_avoiding_line_search(self):
         def fun(v):
-            return np.random.rand(1)
+            return np.random.rand(1)[0]
 
         energy, vr = igl.flip_avoiding_line_search(
             self.f1, self.v1[:, :2], self.v1[:, :2], fun, 10.0)
@@ -2484,20 +2486,22 @@ def test_flip_edge(self):
 
     def test_AABB(self):
         tree = igl.AABB_f64_3()
-        tree.init(self.v1,self.f1)
-        bc = igl.barycenter(self.v1,self.f1)
-        sqrD = tree.squared_distance(self.v1,self.f1,bc)
+        v1_f = np.asarray(self.v1, order='F')
+        f1_f = np.asarray(self.f1, order='F')
+        tree.init(v1_f,f1_f)
+        bc = igl.barycenter(v1_f,f1_f)
+        sqrD = tree.squared_distance(v1_f,f1_f,bc)
         self.assertTrue(sqrD.shape[0] == bc.shape[0])
         self.assertTrue(np.max(sqrD) <= 1e-16)
-        sqrD,I,C = tree.squared_distance(self.v1,self.f1,bc,return_index=True,return_closest_point=True)
+        sqrD,I,C = tree.squared_distance(v1_f,f1_f,bc,return_index=True,return_closest_point=True)
         self.assertTrue(sqrD.shape[0] == bc.shape[0])
         self.assertTrue(I.shape[0] == bc.shape[0])
         self.assertTrue(C.shape == bc.shape)
 
     def test_in_element_3(self):
-        V = np.array([ [0.,0,0], [1,0,0], [0,1,0], [0,0,1], [1,1,1]],dtype='float64')
-        T = np.array([[0,1,2,3],[4,3,2,1]],dtype='int32')
-        Q = np.array([[0.1,0.1,0.1],[0.9,0.9,0.9]],dtype='float64')
+        V = np.array([ [0.,0,0], [1,0,0], [0,1,0], [0,0,1], [1,1,1]],dtype='float64',order='f')
+        T = np.array([[0,1,2,3],[4,3,2,1]],dtype='int32',order='f')
+        Q = np.array([[0.1,0.1,0.1],[0.9,0.9,0.9]],dtype='float64',order='f')
         tree = igl.AABB_f64_3()
         tree.init(V,T)
         I = igl.in_element_3(V,T,Q,tree)
@@ -2506,17 +2510,16 @@ def test_in_element_3(self):
         self.assertTrue(I[1] == 1)
 
     def test_in_element_2(self):
-        V = np.array([ [0.,0], [1,0], [0,1], [1,1]],dtype='float64')
-        F = np.array([[0,1,2],[2,1,3]],'int32')
-        Q = np.array([[0.1,0.1],[0.9,0.9]],dtype='float64')
+        V = np.array([ [0.,0], [1,0], [0,1], [1,1]],dtype='float64',order='f')
+        F = np.array([[0,1,2],[2,1,3]],'int32',order='f')
+        Q = np.array([[0.1,0.1],[0.9,0.9]],dtype='float64',order='f')
         tree = igl.AABB_f64_2()
         tree.init(V,F)
         I = igl.in_element_2(V,F,Q,tree)
         self.assertTrue(I.shape[0] == Q.shape[0])
         self.assertTrue(I[0] == 0)
         self.assertTrue(I[1] == 1)
 
-
     def test_triangulate(self):
         V = np.array([[0,0],[1,0],[1,1],[0,1]],dtype='float64')
         E = np.array([[0,1],[1,2],[2,3],[3,0]])