Merge pull request libigl#275 from libigl/alecjacobson/lipschitz-octree [skip ci]

lipschitz_octree and unique_sparse_voxel_corners

Author: alecjacobson

CMakeLists.txt

@@ -45,7 +45,7 @@ option(LIBIGL_RESTRICTED_TRIANGLE  "Build target igl_restricted::triangle" ON)
 FetchContent_Declare(
   libigl
   GIT_REPOSITORY https://github.com/libigl/libigl.git
-  GIT_TAG cf9ed7f492209590c42dc7247281dfdfb6618487
+  GIT_TAG 678e1fff76815e0c4c5d1f025ee2129181cc7d86
 )
 FetchContent_MakeAvailable(libigl)
 

pyproject.toml

@@ -14,7 +14,7 @@ build-backend = "scikit_build_core.build"
 
 [project]
 name = "libigl"
-version = "2.6.1.dev0"
+version = "2.6.2.dev0"
 description = "libigl: A simple C++ geometry processing library"
 readme = "README.md"
 requires-python = ">=3.8"

src/lipschitz_octree.cpp

@@ -0,0 +1,62 @@
+#include "default_types.h"
+#include <igl/lipschitz_octree.h>
+#include <nanobind/nanobind.h>
+#include <nanobind/ndarray.h>
+#include <nanobind/eigen/dense.h>
+#include <nanobind/stl/tuple.h>
+#include <nanobind/stl/function.h>
+
+namespace nb = nanobind;
+using namespace nb::literals;
+
+namespace pyigl
+{
+  // Hopefully this will avoid a copy by matching the type in lipschitz_octree_prune
+  using MatrixNX3R = Eigen::Matrix<Numeric, Eigen::Dynamic, 3, Eigen::RowMajor>;
+  auto lipschitz_octree(
+    const nb::DRef<const Eigen::VectorXN> &origin,
+    const Numeric h0,
+    const Integer max_depth,
+    const nb::typed<
+      nb::callable,
+      Eigen::VectorXN(const nb::DRef<const MatrixNX3R> &)> & udf)
+  {
+    Eigen::MatrixXI ijk;
+    // Instead of MatrixXN if use use MatrixNX3R we might avoid a copy
+    const std::function<Eigen::VectorXN(const MatrixNX3R &)> udf_wrapper =
+      [&](const MatrixNX3R &Q) -> Eigen::VectorXN
+    {
+      return nb::cast<Eigen::VectorXN>(udf(Q));
+    };
+    igl::lipschitz_octree<true>(origin, h0, max_depth, udf_wrapper, ijk);
+    return ijk;
+  }
+}
+
+// Bind the wrapper to the Python module
+void bind_lipschitz_octree(nb::module_ &m)
+{
+  m.def(
+    "lipschitz_octree",
+    &pyigl::lipschitz_octree, 
+    "origin"_a, 
+    "h0"_a,
+    "max_depth"_a,
+    "udf"_a,
+R"(Given a minimum corner position (origin) and a side length (h0) and a
+maximum depth (max_depth), determine the possible active leaf octree cells
+based on an one-Lipschitz non-negative function to a level set (e.g.,
+"unsigned distance function").
+
+  @param[in] origin  3-vector of root cell origin (minimum corner)
+  @param[in] h0   side length of root cell
+  @param[in] max_depth  maximum depth of octree (root is depth=0)
+  @param[in] udf  1-Lipschitz function of (unsigned) distance to level set to a
+    list of batched query points
+  @param[out] ijk #ijk by 3 list of octree leaf cell minimum corner
+    subscripts
+)");
+}
+
+
+

src/marching_cubes.cpp

@@ -77,7 +77,7 @@ EV = np.array([[k & 0xFFFFFFFF, k >> 32, v] for k, v in E2V.items()], dtype=np.i
     "S"_a, 
     "GV"_a,
     "GI"_a,
-    "isovalue"_a=0,
+    "isovalue"_a=0.0,
 R"(Performs marching cubes reconstruction on a grid defined by values, and
 points, and generates a mesh defined by vertices and faces
 

src/offset_surface.cpp

@@ -10,23 +10,25 @@ namespace nb = nanobind;
 using namespace nb::literals;
 
 namespace pyigl {
-auto offset_surface(const nb::DRef<const Eigen::MatrixXN> &V,
-                    const nb::DRef<const Eigen::MatrixXI> &F,
-                    const Numeric isolevel, const Integer s,
-                    const igl::SignedDistanceType signed_distance_type) {
-
-  Eigen::MatrixXN SV;
-  Eigen::MatrixXI SF;
-  Eigen::MatrixXN GV;
-
-  Eigen::VectorXI side;
-  Eigen::MatrixXN so;
-
-  igl::offset_surface(V, F, isolevel, s, signed_distance_type, SV, SF, GV, side,
-                      so);
-
-  return std::make_tuple(SV, SF, GV, side, so);
-}
+  auto offset_surface(
+    const nb::DRef<const Eigen::MatrixXN> &V,
+    const nb::DRef<const Eigen::MatrixXI> &F,
+    const Numeric isolevel, 
+    const Integer s,
+    const igl::SignedDistanceType signed_distance_type) 
+  {
+    Eigen::MatrixXN SV;
+    Eigen::MatrixXI SF;
+    Eigen::MatrixXN GV;
+  
+    Eigen::VectorXI side;
+    Eigen::MatrixXN so;
+    igl::offset_surface(
+      V, F, isolevel, s, signed_distance_type, SV, SF, GV, side,
+      so);
+  
+    return std::make_tuple(SV, SF, GV, side, so);
+  }
 } // namespace pyigl
 
 // Bind the wrapper to the Python module

src/unique_sparse_voxel_corners.cpp

@@ -0,0 +1,60 @@
+#include "default_types.h"
+#include <igl/unique_sparse_voxel_corners.h>
+#include <nanobind/nanobind.h>
+#include <nanobind/ndarray.h>
+#include <nanobind/eigen/dense.h>
+#include <nanobind/stl/tuple.h>
+
+namespace nb = nanobind;
+using namespace nb::literals;
+
+namespace pyigl
+{
+  auto unique_sparse_voxel_corners(
+    const nb::DRef<const Eigen::VectorXN> &origin,
+    const Numeric h0,
+    const Integer depth,
+    const nb::DRef<const Eigen::MatrixXI> &ijk)
+  {
+    Eigen::MatrixXI unique_ijk;
+    Eigen::MatrixXI J;
+    Eigen::MatrixXN unique_corners;
+    
+    igl::unique_sparse_voxel_corners(
+      origin, h0, depth, ijk, unique_ijk, J, unique_corners);
+    
+    return std::make_tuple(unique_ijk, J, unique_corners);
+  }
+}
+
+// Bind the wrapper to the Python module
+void bind_unique_sparse_voxel_corners(nb::module_ &m)
+{
+  m.def(
+    "unique_sparse_voxel_corners",
+    &pyigl::unique_sparse_voxel_corners, 
+    "origin"_a, 
+    "h0"_a,
+    "depth"_a,
+    "ijk"_a,
+    R"(
+Give a list of octree cells subscripts (ijk) (minimum corners) at a given depth,
+determine a unique list of subscripts to all incident corners of those
+cells (de-replicating shared corners).
+
+ @param[in] origin  3-vector of root cell minimum
+ @param[in] h0   side length of current depth level
+ @param[in] depth  current depth (single root cell is depth = 0)
+ @param[in] ijk #ijk by 3 list of octree leaf cell minimum corner
+   subscripts
+ @param[out] unique_ijk #unique_ijk by 3 list of unique corner subscripts
+ @param[out] J  #ijk by 8 list of indices into unique_ijk in yxz binary
+   counting order
+ @param[out] unique_corners #unique_ijk by 3 list of unique corner
+   positions
+    )");
+
+}
+
+
+

tests/test_all.py

@@ -575,6 +575,18 @@ def test_misc():
     R = igl.oriented_bounding_box(V)
     R = igl.oriented_bounding_box(V,n=100,minimize_type=igl.ORIENTED_BOUNDING_BOX_MINIMIZE_SURFACE_AREA)
 
-
-
-
+def test_octree():
+
+    def sdf_sphere(Q):
+        return np.linalg.norm(Q,axis=1) - 0.5
+    def udf_sphere(Q):
+        return np.abs(sdf_sphere(Q))
+
+    origin = np.array([-1,-1,-1],dtype=np.float64)
+    h0 = 2.0
+    max_depth = 4
+    ijk = igl.lipschitz_octree(origin,h0,max_depth,udf_sphere)
+    h = h0 / (2**max_depth)
+    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)