Merge pull request #182 from libigl/alecjacobson/AABB

Basic AABB support; in_element; (req C++17)

Author: alecjacobson

.github/workflows/wheels.yml

@@ -57,7 +57,7 @@ jobs:
       CIBW_ENVIRONMENT: "MAX_JOBS=${{ matrix.os.runs-on == 'macos-latest' && 3 || 2 }}"
       # Why universal2 here? It's not included above in CIBW_BUILD
       CIBW_ARCHS_MACOS: "x86_64 arm64 universal2"
-      CIBW_ENVIRONMENT_MACOS: "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' || '' }}\""
+      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' || '' }}\""
 
 
     steps:

CMakeLists.txt

@@ -93,7 +93,10 @@ pyigl_include("copyleft" "cgal")
 pyigl_include("restricted" "triangle")
 
 
-add_library(pyigl_classes MODULE classes/classes.cpp)
+file(GLOB PYIGL_CLASSES_SOURCES classes/*.cpp)
+add_library(pyigl_classes MODULE ${PYIGL_CLASSES_SOURCES})
+# std::variant
+target_compile_features(pyigl_classes  PRIVATE cxx_std_17)
 target_link_libraries(pyigl_classes PRIVATE npe igl::core)
 target_link_libraries(pyigl_classes PRIVATE pybind11::module)
 set_target_properties(pyigl_classes PROPERTIES PREFIX "${PYTHON_MODULE_PREFIX}" SUFFIX "${PYTHON_MODULE_EXTENSION}")

classes/AABB.cpp

@@ -0,0 +1,62 @@
+#include <igl/AABB.h>
+#include <npe.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <variant>
+#include <variant>
+
+
+namespace py = pybind11;
+template <int DIM>
+void init_AABB(py::module_ &m)
+{
+  using AABB_f64_DIM = igl::AABB<Eigen::MatrixXd,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)
+      {
+        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);
+        }
+      },
+      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_ &);
+template void init_AABB<3>(py::module_ &);

classes/classes.cpp

@@ -9,8 +9,15 @@
 
 namespace py = pybind11;
 
+// Forward declaration
+template <int DIM>
+void init_AABB(py::module_ &);
+
 PYBIND11_MODULE(pyigl_classes, m)
 {
+  init_AABB<2>(m);
+  init_AABB<3>(m);
+
   py::class_<igl::ARAPData>(m, "ARAP")
       .def(py::init([](Eigen::MatrixXd &v, Eigen::MatrixXi &f, int dim, Eigen::MatrixXi &b,
                        const int energy_type, const bool with_dynamics, const double h, const double ym, const int max_iter) {

src/in_element.cpp

@@ -0,0 +1,54 @@
+#include <common.h>
+#include <npe.h>
+#include <typedefs.h>
+#include <igl/in_element.h>
+#include <igl/AABB.h>
+
+using AABB_f64_3 = igl::AABB<Eigen::MatrixXd,3>;
+const char* ds_in_element = R"igl_Qu8mg5v7(
+Determine whether each point in a list of points is in the elements of a mesh.
+ 
+Parameters:
+------
+V : #V by dim list of mesh vertex positions.
+Ele : #Ele by dim+1 list of mesh indices into #V.
+Q  : #Q by dim list of query point positions
+aabb :  axis-aligned bounding box tree object (see AABB.h)
+
+Returns:
+-------
+I : #Q list of indices into Ele of first containing element (-1 means no
+containing element)
+)igl_Qu8mg5v7";
+npe_function(in_element_3)
+npe_doc(  ds_in_element)
+npe_arg(V, Eigen::MatrixXd)
+npe_arg(Ele, Eigen::MatrixXi)
+npe_arg(Q, Eigen::MatrixXd)
+npe_arg(aabb, AABB_f64_3)
+npe_begin_code()
+
+  Eigen::VectorXi I;
+  igl::in_element(V,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)
+npe_arg(Ele, Eigen::MatrixXi)
+npe_arg(Q, Eigen::MatrixXd)
+npe_arg(aabb, AABB_f64_2)
+npe_begin_code()
+
+  Eigen::VectorXi I;
+  igl::in_element(V,Ele,Q,aabb,I);
+  return npe::move(I);
+
+npe_end_code()
+
+
+

tests/test_basic.py

@@ -2470,6 +2470,41 @@ def test_flip_edge(self):
                         emap.dtype == self.f1.dtype)
         self.assertTrue(np.array(ue2e).dtype == self.f1.dtype)
 
+    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)
+        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)
+        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')
+        tree = igl.AABB_f64_3()
+        tree.init(V,T)
+        I = igl.in_element_3(V,T,Q,tree)
+        self.assertTrue(I.shape[0] == Q.shape[0])
+        self.assertTrue(I[0] == 0)
+        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')
+        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]])
@@ -2489,6 +2524,7 @@ def test_triangulate(self):
         self.assertTrue(E2.shape == E.shape)
         self.assertTrue(EM2.shape == EM.shape)
 
+
     # copyleft.cgal
     def test_convex_hull(self):
         V = np.array([[0,0,0],[1,0,0],[0,1,0],[0,0,1]],dtype="float64")