Lipschitz octree pruning (#2488) [skip ci]

* lipshitz_octree tutorial example

* better comments

* revert to shallow default

* better names and documentation, dont return I

* long long -> std::int64_t

* asdd tutorail chapter 10

* precious windows

* int64_t -> std::int64_t

* batch version

* fix template chain from overzealous index typing

* bump to c++17 (if constexpr)

.github/workflows/continuous.yml

@@ -114,6 +114,7 @@ jobs:
         {static: OFF, tutorials:  ON, tests: OFF, selected_tutorial: 7},
         {static: OFF, tutorials:  ON, tests: OFF, selected_tutorial: 8},
         {static: OFF, tutorials:  ON, tests: OFF, selected_tutorial: 9},
+        {static: OFF, tutorials:  ON, tests: OFF, selected_tutorial: 10},
       ]
     steps:
       - name: Checkout repository
@@ -166,6 +167,7 @@ jobs:
             -DLIBIGL_TUTORIALS_CHAPTER7=${{ (matrix.build-params.selected_tutorial == 'NONE' || matrix.build-params.selected_tutorial == '7') && 'ON' || 'OFF' }} ^
             -DLIBIGL_TUTORIALS_CHAPTER8=${{ (matrix.build-params.selected_tutorial == 'NONE' || matrix.build-params.selected_tutorial == '8') && 'ON' || 'OFF' }} ^
             -DLIBIGL_TUTORIALS_CHAPTER9=${{ (matrix.build-params.selected_tutorial == 'NONE' || matrix.build-params.selected_tutorial == '9') && 'ON' || 'OFF' }} ^
+            -DLIBIGL_TUTORIALS_CHAPTER10=${{ (matrix.build-params.selected_tutorial == 'NONE' || matrix.build-params.selected_tutorial == '10') && 'ON' || 'OFF' }} ^
             -B build ^
             -S .
           cmake --build build -j2

cmake/igl/igl_add_library.cmake

@@ -33,8 +33,8 @@ function(igl_add_library module_name)
         target_compile_definitions(${module_name} ${IGL_SCOPE} -DIGL_STATIC_LIBRARY)
     endif()
 
-    # C++11 features
-    target_compile_features(${module_name} ${IGL_SCOPE} cxx_std_11)
+    # C++17 features
+    target_compile_features(${module_name} ${IGL_SCOPE} cxx_std_17)
 
     if(LIBIGL_WARNINGS_AS_ERRORS)
       target_compile_options(${module_name} PRIVATE -Wall -Wextra -Wpedantic -Wno-sign-compare -Werror -Wno-gnu -Wno-unknown-pragmas)

include/igl/is_vertex_manifold.cpp

@@ -22,7 +22,7 @@ IGL_INLINE bool igl::is_vertex_manifold(
 {
   assert(F.cols() == 3 && "F must contain triangles");
   typedef typename DerivedF::Scalar Index;
-  typedef typename DerivedF::Index FIndex;
+  using FIndex = int;
   const Index n = F.maxCoeff()+1;
   std::vector<std::vector<std::vector<FIndex > > > TT;
   std::vector<std::vector<std::vector<FIndex > > > TTi;

include/igl/lipschitz_octree.cpp

@@ -0,0 +1,97 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2025 Alec Jacobson <[email protected]>
+//
+// 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/.
+#include "lipschitz_octree.h"
+#include "lipschitz_octree_prune.h"
+#include "find.h"
+#include "matlab_format.h"
+#include <cassert>
+#include <iostream>
+
+template <
+  bool batched,
+  typename Derivedorigin,
+  typename Func,
+  typename Derivedijk
+    >
+IGL_INLINE void igl::lipschitz_octree(
+  const Eigen::MatrixBase<Derivedorigin> & origin,
+  const typename Derivedorigin::Scalar h0,
+  const int max_depth,
+  const Func & udf,
+  Eigen::PlainObjectBase<Derivedijk> & ijk_out)
+{
+  using Scalar = typename Derivedorigin::Scalar;
+  using RowVectorS3 = Eigen::Matrix<Scalar,1,3>;
+  using MatrixSX8R = Eigen::Matrix<Scalar,Eigen::Dynamic,8,Eigen::RowMajor>;
+  using MatrixiX3R = Eigen::Matrix<int,Eigen::Dynamic,3,Eigen::RowMajor>;
+
+  // static assert to ensure that Derivedorigin is a vector and the
+  // non-singleton dimension is 3 or Eigen::Dynamic
+  static_assert(
+    (Derivedorigin::RowsAtCompileTime == 1 && (
+      Derivedorigin::ColsAtCompileTime == 3 ||
+      Derivedorigin::ColsAtCompileTime == Eigen::Dynamic)) ||
+    (Derivedorigin::ColsAtCompileTime == 1 && (
+      Derivedorigin::RowsAtCompileTime == 3 ||
+      Derivedorigin::RowsAtCompileTime == Eigen::Dynamic)),
+    "Derivedorigin must be a vector with 3 or Eigen::Dynamic dimensions");
+  // dynamic assert that the origin is a 3D vector
+  assert((origin.rows() == 3 || origin.cols() == 3) && origin.size() == 3 &&
+    "origin must be a 3D vector");
+
+  MatrixiX3R ijk(1,3);
+  ijk<<0,0,0;
+
+  for(int depth = 0;depth<=max_depth;depth++)
+  {
+    if(ijk.rows() == 0)
+    {
+      // no more cells to refine
+      break;
+    }
+    const Scalar h = h0 / (1 << depth);
+    MatrixiX3R ijk_next;
+    MatrixiX3R ijk_maybe;
+    igl::lipschitz_octree_prune<batched>(origin,h0,depth,udf,ijk,ijk_maybe);
+    if(depth == max_depth)
+    {
+      // sad copy
+      ijk_out = ijk_maybe.template cast<typename Derivedijk::Scalar>();
+      return;
+    }else
+    {
+      const MatrixiX3R ijk_maybe_2 = ijk_maybe * 2;
+
+      ijk_next.resize(ijk_maybe_2.rows()*8,3);
+      for(int c = 0;c<ijk_maybe_2.rows();c++)
+      {
+        for(int i = 0;i<8;i++)
+        {
+          const int k = i;
+          ijk_next.row(c*8+k) =
+            ijk_maybe_2.row(c) +
+            // This order shouldn't really matter, though it seems nice if it
+            // matches above.
+            Eigen::RowVector3i((i&2) ? 1 : 0, (i&4) ? 1 : 0, (i&1) ? 1 : 0);
+        }
+      }
+    }
+    if(depth == max_depth)
+    {
+      break;
+    }
+    ijk = ijk_next;
+  }
+  assert(false && "Should never reach here");
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::lipschitz_octree<false,Eigen::Matrix<double, 1, 3, 1, 1, 3>, std::function<double (Eigen::Matrix<double, 1, 3, 1, 1, 3> const&)>, Eigen::Matrix<int, -1, 3, 1, -1, 3>>(Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3>> const&, Eigen::Matrix<double, 1, 3, 1, 1, 3>::Scalar, int, std::function<double (Eigen::Matrix<double, 1, 3, 1, 1, 3> const&)> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3>>&);
+template void igl::lipschitz_octree<true, Eigen::Matrix<double, 1, 3, 1, 1, 3>, std::function<Eigen::Matrix<double, -1, 1, 0, -1, 1> (Eigen::Matrix<double, -1, 3, 1, -1, 3> const&)>, Eigen::Matrix<int, -1, 3, 1, -1, 3>>(Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3>> const&, Eigen::Matrix<double, 1, 3, 1, 1, 3>::Scalar, int, std::function<Eigen::Matrix<double, -1, 1, 0, -1, 1> (Eigen::Matrix<double, -1, 3, 1, -1, 3> const&)> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3>>&);
+#endif

include/igl/lipschitz_octree.h

@@ -0,0 +1,46 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2025 Alec Jacobson <[email protected]>
+//
+// 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/.
+#ifndef IGL_LIPSCHITZ_OCTREE_H
+#define IGL_LIPSCHITZ_OCTREE_H
+
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl 
+{
+  /// 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
+  ///   @param[out] ijk #ijk by 3 list of octree leaf cell minimum corner
+  ///     subscripts
+  template <
+    bool batched=false,
+    typename Derivedorigin,
+    typename Func,
+    typename Derivedijk
+      >
+  IGL_INLINE void lipschitz_octree(
+    const Eigen::MatrixBase<Derivedorigin> & origin,
+    const typename Derivedorigin::Scalar h0,
+    const int max_depth,
+    const Func & udf,
+    Eigen::PlainObjectBase<Derivedijk> & ijk);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#    include "lipschitz_octree.cpp"
+#endif
+
+#endif
+

include/igl/lipschitz_octree_prune.cpp

@@ -0,0 +1,109 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2025 Alec Jacobson <[email protected]>
+//
+// 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/.
+#include "lipschitz_octree_prune.h"
+#include "unique_sparse_voxel_corners.h"
+#include "find.h"
+#include "matlab_format.h"
+#include "parallel_for.h"
+#include <cassert>
+#include <iostream>
+
+
+template <
+  bool batched,
+  typename Derivedorigin,
+  typename Func,
+  typename Derivedijk,
+  typename Derivedijk_maybe
+    >
+IGL_INLINE void igl::lipschitz_octree_prune(
+  const Eigen::MatrixBase<Derivedorigin> & origin,
+  const typename Derivedorigin::Scalar h0,
+  const int depth,
+  const Func & udf,
+  const Eigen::MatrixBase<Derivedijk> & ijk,
+  Eigen::PlainObjectBase<Derivedijk_maybe> & ijk_maybe)
+{
+  // static assert to ensure that Derivedorigin is a vector and the
+  // non-singleton dimension is 3 or Eigen::Dynamic
+  static_assert(
+    (Derivedorigin::RowsAtCompileTime == 1 && (
+      Derivedorigin::ColsAtCompileTime == 3 ||
+      Derivedorigin::ColsAtCompileTime == Eigen::Dynamic)) ||
+    (Derivedorigin::ColsAtCompileTime == 1 && (
+      Derivedorigin::RowsAtCompileTime == 3 ||
+      Derivedorigin::RowsAtCompileTime == Eigen::Dynamic)),
+    "Derivedorigin must be a vector with 3 or Eigen::Dynamic dimensions");
+  // dynamic assert that the origin is a 3D vector
+  assert((origin.rows() == 3 || origin.cols() == 3) && origin.size() == 3 &&
+    "origin must be a 3D vector");
+
+  using Scalar = typename Derivedorigin::Scalar;
+  using RowVectorS3 = Eigen::Matrix<Scalar,1,3>;
+  using MatrixSX8R = Eigen::Matrix<typename Derivedorigin::Scalar,Eigen::Dynamic,8,Eigen::RowMajor>;
+  using MatrixSX3R = Eigen::Matrix<typename Derivedorigin::Scalar,Eigen::Dynamic,3,Eigen::RowMajor>;
+  using MatrixiX3R = Eigen::Matrix<int,Eigen::Dynamic,3,Eigen::RowMajor>;
+  // h0 is already the h at this depth.
+  const Scalar h = h0 / (1 << depth);
+
+  Eigen::Matrix<int,Eigen::Dynamic,8,Eigen::RowMajor> J;
+  MatrixiX3R unique_ijk;
+  MatrixSX3R unique_corner_positions;
+  igl::unique_sparse_voxel_corners(origin,h0,depth,ijk,unique_ijk,J,unique_corner_positions);
+
+  // Effectively a batched call to udf
+  Eigen::Array<bool,Eigen::Dynamic,1> big(unique_corner_positions.rows()); 
+  // Requires C++17
+  if constexpr (batched)
+  {
+    Eigen::Matrix<Scalar,Eigen::Dynamic,1> u = udf(unique_corner_positions);
+    for(int i = 0;i<unique_corner_positions.rows();i++)
+    {
+      assert(u(i) >= 0 && "udf must be non-negative for lipschitz_octree_prune");
+      big(i) = (u(i) > h * std::sqrt(3));
+    }
+  }else
+  {
+    //for(int u = 0;u<unique_corner_positions.rows();u++)
+    igl::parallel_for(
+      unique_corner_positions.rows(),
+      [&](const int i)
+      {
+        // evaluate the function at the corner
+        const RowVectorS3 corner = unique_corner_positions.row(i);
+        const Scalar u = udf(corner);
+        assert(u >= 0 && "udf must be non-negative for lipschitz_octree_prune");
+        big(i) = (u > h * std::sqrt(3));
+      },
+      1000);
+  }
+
+  ijk_maybe.resize(ijk.rows(),3);
+  int k = 0;
+  for(int c = 0;c<ijk.rows();c++)
+  {
+    bool empty = false;
+    for(int i = 0;i<8;i++)
+    {
+      empty = empty || big(J(c,i));
+    }
+    bool maybe = !empty;
+    if(maybe)
+    {
+      ijk_maybe.row(k++) = ijk.row(c);
+    }
+  }
+  ijk_maybe.conservativeResize(k,3);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::lipschitz_octree_prune<false,Eigen::Matrix<double, 1, 3, 1, 1, 3>, std::function<double (Eigen::Matrix<double, 1, 3, 1, 1, 3> const&)>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>>(Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3>> const&, Eigen::Matrix<double, 1, 3, 1, 1, 3>::Scalar, int, std::function<double (Eigen::Matrix<double, 1, 3, 1, 1, 3> const&)> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3>> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3>>&);
+template void igl::lipschitz_octree_prune<true, Eigen::Matrix<double, 1, 3, 1, 1, 3>, std::function<Eigen::Matrix<double, -1, 1, 0, -1, 1> (Eigen::Matrix<double, -1, 3, 1, -1, 3> const&)>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>>(Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3>> const&, Eigen::Matrix<double, 1, 3, 1, 1, 3>::Scalar, int, std::function<Eigen::Matrix<double, -1, 1, 0, -1, 1> (Eigen::Matrix<double, -1, 3, 1, -1, 3> const&)> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3>> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3>>&);
+#endif
+

include/igl/lipschitz_octree_prune.h

@@ -0,0 +1,52 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2025 Alec Jacobson <[email protected]>
+//
+// 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/.
+#ifndef IGL_LIPSCHITZ_OCTREE_PRUNE_H
+#define IGL_LIPSCHITZ_OCTREE_PRUNE_H
+
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl 
+{
+  /// Given a 1-Lipschitz non-negative function to a level set (e.g., "unsigned
+  /// distance function") and a set of query octree cells at a current depth
+  /// given by the minimum corners subscripts in ijk (where root is (0,0,0)),
+  /// determine which of these cells maybe still contains the level set based on
+  /// checking a simple sufficient culling condition at the corners.
+  ///
+  ///   @param[in] origin  3-vector of root cell origin (minimum corner)
+  ///   @param[in] h0   side length of root cell
+  ///   @param[in] depth  current depth (single root cell is depth=0)
+  ///   @param[in] udf  1-Lipschitz function of (unsigned) distance to level set
+  ///   @param[in] ijk #ijk by 3 list of octree leaf cell minimum corner
+  ///   subscripts
+  ///   @param[out] ijk_maybe #ijk_maybe by 3 list of octree leaf cell indices
+  template <
+    bool batched=false,
+    typename Derivedorigin,
+    typename Func,
+    typename Derivedijk,
+    typename Derivedijk_maybe
+      >
+  IGL_INLINE void lipschitz_octree_prune(
+    const Eigen::MatrixBase<Derivedorigin> & origin,
+    const typename Derivedorigin::Scalar h0,
+    const int depth,
+    const Func & udf,
+    const Eigen::MatrixBase<Derivedijk> & ijk,
+    Eigen::PlainObjectBase<Derivedijk_maybe> & ijk_maybe);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#    include "lipschitz_octree_prune.cpp"
+#endif
+
+#endif
+
+
+

include/igl/list_to_matrix.cpp

@@ -155,6 +155,7 @@ IGL_INLINE bool igl::list_to_matrix(const std::vector<T > & V,Eigen::PlainObject
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+template bool igl::list_to_matrix<std::int64_t, Eigen::Matrix<std::int64_t, -1, 1, 0, -1, 1>>(std::vector<std::int64_t, std::allocator<std::int64_t>> const&, Eigen::PlainObjectBase<Eigen::Matrix<std::int64_t, -1, 1, 0, -1, 1>>&);
 // generated by autoexplicit.sh
 template bool igl::list_to_matrix<int, Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> >(std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> >&);
 // generated by autoexplicit.sh

include/igl/march_cube.cpp

@@ -133,11 +133,14 @@ IGL_INLINE void igl::march_cube(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::march_cube<Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>>, double, long, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>> const&, Eigen::Matrix<double, 8, 1, 0, 8, 1> const&, Eigen::Matrix<long, 8, 1, 0, 8, 1> const&, double const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3>>&, long&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3>>&, long&, std::unordered_map<std::int64_t, int, std::hash<std::int64_t>, std::equal_to<std::int64_t>, std::allocator<std::pair<std::int64_t const, int>>>&);
 template void igl::march_cube<Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> >, float, unsigned int, Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<float, 8, 1, 0, 8, 1> const&, Eigen::Matrix<unsigned int, 8, 1, 0, 8, 1> const&, float const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&, unsigned int&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> >&, unsigned int&, std::unordered_map<std::int64_t, int, std::hash<std::int64_t>, std::equal_to<std::int64_t>, std::allocator<std::pair<std::int64_t const, int> > >&);
 template void igl::march_cube<Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >, double, unsigned int, Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<double, 8, 1, 0, 8, 1> const&, Eigen::Matrix<unsigned int, 8, 1, 0, 8, 1> const&, double const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&, unsigned int&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> >&, unsigned int&, std::unordered_map<std::int64_t, int, std::hash<std::int64_t>, std::equal_to<std::int64_t>, std::allocator<std::pair<std::int64_t const, int> > >&);
 template void igl::march_cube<Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >, double, long, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<double, 8, 1, 0, 8, 1> const&, Eigen::Matrix<long, 8, 1, 0, 8, 1> const&, double const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, long&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, long&, std::unordered_map<std::int64_t, int, std::hash<std::int64_t>, std::equal_to<std::int64_t>, std::allocator<std::pair<std::int64_t const, int> > >&);
 template void igl::march_cube<Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >, double, unsigned int, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<double, 8, 1, 0, 8, 1> const&, Eigen::Matrix<unsigned int, 8, 1, 0, 8, 1> const&, double const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, unsigned int&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, unsigned int&, std::unordered_map<std::int64_t, int, std::hash<std::int64_t>, std::equal_to<std::int64_t>, std::allocator<std::pair<std::int64_t const, int> > >&);
 #ifdef WIN32
+template void __cdecl igl::march_cube<class Eigen::MatrixBase<class Eigen::Matrix<double,-1,3,1,-1,3> >,double,__int64,class Eigen::Matrix<double,-1,3,0,-1,3>,class Eigen::Matrix<int,-1,3,0,-1,3> >(class Eigen::MatrixBase<class Eigen::Matrix<double,-1,3,1,-1,3> > const &,class Eigen::Matrix<double,8,1,0,8,1> const &,class Eigen::Matrix<__int64,8,1,0,8,1> const &,double const &,class Eigen::PlainObjectBase<class Eigen::Matrix<double,-1,3,0,-1,3> > &,__int64 &,class Eigen::PlainObjectBase<class Eigen::Matrix<int,-1,3,0,-1,3> > &,__int64 &,class std::unordered_map<__int64,int,struct std::hash<__int64>,struct std::equal_to<__int64>,class std::allocator<struct std::pair<__int64 const ,int> > > &);
 template void __cdecl igl::march_cube<class Eigen::MatrixBase<class Eigen::Matrix<double,-1,-1,0,-1,-1> >,double,__int64,class Eigen::Matrix<double,-1,-1,0,-1,-1>,class Eigen::Matrix<int,-1,-1,0,-1,-1> >(class Eigen::MatrixBase<class Eigen::Matrix<double,-1,-1,0,-1,-1> > const &,class Eigen::Matrix<double,8,1,0,8,1> const &,class Eigen::Matrix<__int64,8,1,0,8,1> const &,double const &,class Eigen::PlainObjectBase<class Eigen::Matrix<double,-1,-1,0,-1,-1> > &,__int64 &,class Eigen::PlainObjectBase<class Eigen::Matrix<int,-1,-1,0,-1,-1> > &,__int64 &,class std::unordered_map<__int64,int,struct std::hash<__int64>,struct std::equal_to<__int64>,class std::allocator<struct std::pair<__int64 const ,int> > > &);
 #endif
 #endif 

include/igl/marching_cubes.cpp

@@ -149,6 +149,7 @@ IGL_INLINE void igl::marching_cubes(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+template void igl::marching_cubes<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 8, 1, -1, 8>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 8, 1, -1, 8>> const&, Eigen::Matrix<double, -1, 1, 0, -1, 1>::Scalar, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3>>&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3>>&);
 // generated by autoexplicit.sh
 template void igl::marching_cubes<Eigen::Matrix<float, -1, 1, 0, -1, 1>, Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, unsigned int, unsigned int, unsigned int, Eigen::Matrix<float, -1, 1, 0, -1, 1>::Scalar, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> >&);
 // generated by autoexplicit.sh

include/igl/marching_cubes.h

@@ -67,7 +67,8 @@ namespace igl
   ///
   /// @param[in] S #S list of scalar field values
   /// @param[in] GV  #S by 3 list of referenced grid vertex positions
-  /// @param[in] GI  #GI by 8 list of grid corner indices into rows of GV
+  /// @param[in] GI  #GI by 8 list of grid corner indices into rows of GV (e.g.,
+  /// as output by igl::sparse_voxel_grid) in y-x-z binary counting order.
   template <
     typename DerivedS, 
     typename DerivedGV, 

include/igl/matrix_to_list.cpp

@@ -53,6 +53,7 @@ IGL_INLINE std::vector<typename DerivedM::Scalar > igl::matrix_to_list(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+template void igl::matrix_to_list<Eigen::Matrix<std::int64_t, -1, 8, 1, -1, 8>>(Eigen::MatrixBase<Eigen::Matrix<std::int64_t, -1, 8, 1, -1, 8>> const&, std::vector<Eigen::Matrix<std::int64_t, -1, 8, 1, -1, 8>::Scalar, std::allocator<Eigen::Matrix<std::int64_t, -1, 8, 1, -1, 8>::Scalar>>&);
 // generated by autoexplicit.sh
 template void igl::matrix_to_list<Eigen::Block<Eigen::Matrix<float, -1, -1, 0, -1, -1>, -1, 1, true> >(Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<float, -1, -1, 0, -1, -1>, -1, 1, true> > const&, std::vector<Eigen::Block<Eigen::Matrix<float, -1, -1, 0, -1, -1>, -1, 1, true>::Scalar, std::allocator<Eigen::Block<Eigen::Matrix<float, -1, -1, 0, -1, -1>, -1, 1, true>::Scalar> >&);
 // generated by autoexplicit.sh

include/igl/signed_distance.h

@@ -264,7 +264,7 @@ namespace igl
   ///     Eigen::VectorXd V, P; //where V is mesh vertices, P are query points
   ///     Eigen::VectorXi F;
   ///     igl::FastWindingNumberBVH fwn_bvh;
-  ///     igl::fast_winding_number(V.cast<float>(), F, 2, fwn_bvh);
+  ///     igl::fast_winding_number(V.cast<float>().eval(), F, 2, fwn_bvh);
   ///     igl::signed_distance_fast_winding_number(P,V,F,tree,fwn_bvh,S);
   ///
   /// @param[in] P  #P by 3 list of query point positions

include/igl/slice.cpp

@@ -207,6 +207,8 @@ IGL_INLINE void igl::slice(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+template void igl::slice<std::int64_t>(std::vector<std::int64_t, std::allocator<std::int64_t>> const&, std::vector<size_t, std::allocator<size_t>> const&, std::vector<std::int64_t, std::allocator<std::int64_t>>&);
+
 template void igl::slice<Eigen::SparseMatrix<bool, 0, int>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::SparseMatrix<bool, 0, int>>(Eigen::SparseMatrix<bool, 0, int> const &, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>> const &, int, Eigen::SparseMatrix<bool, 0, int> &);
 template void igl::slice<Eigen::SparseMatrix<double, 0, int>, Eigen::Array<int, -1, 1, 0, -1, 1>, Eigen::SparseMatrix<double, 0, int> >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::DenseBase<Eigen::Array<int, -1, 1, 0, -1, 1> > const&, int, Eigen::SparseMatrix<double, 0, int>&);
 template void igl::slice<Eigen::SparseMatrix<double, 0, int>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::SparseMatrix<double, 0, int>>(Eigen::SparseMatrix<double, 0, int> const &, Eigen::DenseBase<Eigen::Matrix<double, -1, 1, 0, -1, 1>> const &, int, Eigen::SparseMatrix<double, 0, int> &);
@@ -217,11 +219,10 @@ template void igl::slice<unsigned int>(class std::vector<unsigned int,class std:
 template void igl::slice<float>(class std::vector<float,class std::allocator<float> > const &,class std::vector<unsigned __int64,class std::allocator<unsigned __int64> > const &,class std::vector<float,class std::allocator<float> > &);
 template void igl::slice<__int64>(class std::vector<__int64,class std::allocator<__int64> > const &,class std::vector<unsigned __int64,class std::allocator<unsigned __int64> > const &,class std::vector<__int64,class std::allocator<__int64> > &);
 #endif
-template void igl::slice<unsigned int>(std::vector<unsigned int, std::allocator<unsigned int> > const&, std::vector<unsigned long, std::allocator<unsigned long> > const&, std::vector<unsigned int, std::allocator<unsigned int> >&);
+template void igl::slice<unsigned int>(std::vector<unsigned int, std::allocator<unsigned int> > const&, std::vector<size_t, std::allocator<size_t> > const&, std::vector<unsigned int, std::allocator<unsigned int> >&);
 template void igl::slice<float>(std::vector<float, std::allocator<float> > const&, std::vector<size_t, std::allocator<size_t> > const&, std::vector<float, std::allocator<float> >&);
 template void igl::slice<double>(std::vector<double, std::allocator<double> > const&, std::vector<size_t, std::allocator<size_t> > const&, std::vector<double, std::allocator<double> >&);
 template void igl::slice<int>(std::vector<int, std::allocator<int> > const&, std::vector<size_t, std::allocator<size_t> > const&, std::vector<int, std::allocator<int> >&);
-template void igl::slice<long>(std::vector<long, std::allocator<long> > const&, std::vector<size_t, std::allocator<size_t> > const&, std::vector<long, std::allocator<long> >&);
 #include "SortableRow.h"
 template void igl::slice<igl::SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> > >(std::vector<igl::SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> >, std::allocator<igl::SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> > > > const&, std::vector<size_t, std::allocator<size_t> > const&, std::vector<igl::SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> >, std::allocator<igl::SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> > > >&);
 template void igl::slice<igl::SortableRow<Eigen::Matrix<double, -1, 1, 0, -1, 1> > >(std::vector<igl::SortableRow<Eigen::Matrix<double, -1, 1, 0, -1, 1> >, std::allocator<igl::SortableRow<Eigen::Matrix<double, -1, 1, 0, -1, 1> > > > const&, std::vector<size_t, std::allocator<size_t> > const&, std::vector<igl::SortableRow<Eigen::Matrix<double, -1, 1, 0, -1, 1> >, std::allocator<igl::SortableRow<Eigen::Matrix<double, -1, 1, 0, -1, 1> > > >&);

include/igl/sort.cpp

@@ -321,6 +321,7 @@ if(!ascending)
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+template void igl::sort<std::int64_t>(std::vector<std::int64_t, std::allocator<std::int64_t>> const&, bool, std::vector<std::int64_t, std::allocator<std::int64_t>>&, std::vector<size_t, std::allocator<size_t>>&);
 // generated by autoexplicit.sh
 template void igl::sort<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::DenseBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 // generated by autoexplicit.sh
@@ -363,7 +364,6 @@ template void igl::sort_new<Eigen::Matrix<int, 1, 6, 1, 1, 6>, Eigen::Matrix<int
 template void igl::sort<Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 2, 0, -1, 2> >(Eigen::DenseBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&);
 template void igl::sort<Eigen::Matrix<double, -1, 4, 0, -1, 4>, Eigen::Matrix<double, -1, 4, 0, -1, 4>, Eigen::Matrix<int, -1, 4, 0, -1, 4> >(Eigen::DenseBase<Eigen::Matrix<double, -1, 4, 0, -1, 4> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 4, 0, -1, 4> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 4, 0, -1, 4> >&);
 template void igl::sort<Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::DenseBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-template void igl::sort<long>(std::vector<long, std::allocator<long> > const&, bool, std::vector<long, std::allocator<long> >&, std::vector<size_t, std::allocator<size_t> >&);
 template void igl::sort<Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::DenseBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::sort<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::DenseBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::sort<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::DenseBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);

include/igl/triangle_triangle_adjacency.cpp

@@ -307,6 +307,8 @@ template <
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::triangle_triangle_adjacency<Eigen::Matrix<int, -1, -1, 0, -1, -1>, int, int>(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>> const&, std::vector<std::vector<std::vector<int, std::allocator<int>>, std::allocator<std::vector<int, std::allocator<int>>>>, std::allocator<std::vector<std::vector<int, std::allocator<int>>, std::allocator<std::vector<int, std::allocator<int>>>>>>&, std::vector<std::vector<std::vector<int, std::allocator<int>>, std::allocator<std::vector<int, std::allocator<int>>>>, std::allocator<std::vector<std::vector<int, std::allocator<int>>, std::allocator<std::vector<int, std::allocator<int>>>>>>&);
+// generated by autoexplicit.sh
 template void igl::triangle_triangle_adjacency<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, int, int>(Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, bool, std::vector<std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >, std::allocator<std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > > >&, std::vector<std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >, std::allocator<std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > > >&);
 // generated by autoexplicit.sh
 template void igl::triangle_triangle_adjacency<Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);

include/igl/unique.cpp

@@ -195,6 +195,7 @@ IGL_INLINE void igl::unique(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+template void igl::unique<Eigen::Matrix<std::int64_t, -1, 8, 1, -1, 8>, Eigen::Matrix<std::int64_t, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>>(Eigen::MatrixBase<Eigen::Matrix<std::int64_t, -1, 8, 1, -1, 8>> const&, Eigen::PlainObjectBase<Eigen::Matrix<std::int64_t, -1, 1, 0, -1, 1>>&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>>&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>>&);
 // generated by autoexplicit.sh
 template void igl::unique<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::unique<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
@@ -210,7 +211,6 @@ template void igl::unique<Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int
 template void igl::unique<double>(std::vector<double, std::allocator<double> > const&, std::vector<double, std::allocator<double> >&);
 template void igl::unique<int>(std::vector<int, std::allocator<int> > const&, std::vector<int, std::allocator<int> >&);
 template void igl::unique<int>(std::vector<int, std::allocator<int> > const&, std::vector<int, std::allocator<int> >&, std::vector<size_t, std::allocator<size_t> >&, std::vector<size_t, std::allocator<size_t> >&);
-template void igl::unique<long>(std::vector<long, std::allocator<long> > const&, std::vector<long, std::allocator<long> >&, std::vector<size_t, std::allocator<size_t> >&, std::vector<size_t, std::allocator<size_t> >&);
 #ifdef WIN32
 template void igl::unique<class Eigen::Matrix<int,-1,1,0,-1,1>,class Eigen::Matrix<int,-1,1,0,-1,1>,class Eigen::Matrix<__int64,-1,1,0,-1,1>,class Eigen::Matrix<__int64,-1,1,0,-1,1> >(class Eigen::MatrixBase<class Eigen::Matrix<int,-1,1,0,-1,1> > const &,class Eigen::PlainObjectBase<class Eigen::Matrix<int,-1,1,0,-1,1> > &,class Eigen::PlainObjectBase<class Eigen::Matrix<__int64,-1,1,0,-1,1> > &,class Eigen::PlainObjectBase<class Eigen::Matrix<__int64,-1,1,0,-1,1> > &);
 template void igl::unique<__int64>(class std::vector<__int64,class std::allocator<__int64> > const &,class std::vector<__int64,class std::allocator<__int64> > &,class std::vector<unsigned __int64,class std::allocator<unsigned __int64> > &,class std::vector<unsigned __int64,class std::allocator<unsigned __int64> > &);

include/igl/unique_sparse_voxel_corners.cpp

@@ -0,0 +1,142 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2025 Alec Jacobson <[email protected]>
+//
+// 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/.
+#include "unique_sparse_voxel_corners.h"
+#include "matlab_format.h"
+#include "unique.h"
+#include <cassert>
+#include <iostream>
+
+namespace
+{
+  // Match the yxz binary counting order in marching_cubes/sparse_voxel_grid
+  const int marching_cubes_reoder[] = {1,0,2,3,5,4,6,7};
+}
+
+template <
+  typename Derivedijk,
+  typename Derivedunique_ijk,
+  typename DerivedJ
+    >
+IGL_INLINE void igl::unique_sparse_voxel_corners(
+  const int depth,
+  const Eigen::MatrixBase<Derivedijk> & ijk,
+  Eigen::PlainObjectBase<Derivedunique_ijk> & unique_ijk,
+  Eigen::PlainObjectBase<DerivedJ> & J)
+{
+  using MatrixiX3R = Eigen::Matrix<int,Eigen::Dynamic,3,Eigen::RowMajor>;
+  // slow, precise hashing. Not sure why I didn't just use unique_rows on #ijk*8
+  // by 3...
+  const Eigen::Matrix<std::int64_t,1,3> coeffs(
+    1,
+    //static_cast<std::int64_t>(std::pow(2,depth) + 1), 
+    (1 << depth) + 1,
+    //static_cast<std::int64_t>(std::pow(2,depth) + 1) * (std::pow(2,depth) + 1));
+    ((1 << depth) +1)*((1 << depth) +1));
+  const auto ijk2code = [&coeffs](const int i, const int j, const int k)->std::int64_t
+  {
+    // code = i*(2.^depth + 1)^0 + j*(2.^depth + 1)^1 + k*(2.^depth + 1)^2;
+    // Probably can just use 2.^(depth+1) instead of (2.^depth + 1) and use bit
+    // shifting?
+    const std::int64_t code =
+      static_cast<std::int64_t>(i) * coeffs[0] +
+      static_cast<std::int64_t>(j) * coeffs[1] +
+      static_cast<std::int64_t>(k) * coeffs[2];
+    return code;
+  };
+  const auto code2ijk = [&coeffs](const std::int64_t code, int & i, int & j, int & k)
+  {
+    k = static_cast<int>(code / coeffs[2]);
+    j = static_cast<int>((code - k * coeffs[2]) / coeffs[1]);
+    i = static_cast<int>(code - j * coeffs[1] - k * coeffs[2]);
+  };
+
+  // Should use parallel_for
+  Eigen::Matrix<std::int64_t,Eigen::Dynamic,8,Eigen::RowMajor> codes(ijk.rows(),8);
+  for(int c = 0;c<ijk.rows();c++)
+  {
+    for(int i = 0;i<8;i++)
+    {
+      Eigen::RowVector3i ijk_c(
+        ijk(c,0) + ((i&2) ? 1 : 0), 
+        ijk(c,1) + ((i&4) ? 1 : 0), 
+        ijk(c,2) + ((i&1) ? 1 : 0));
+      const std::int64_t code = ijk2code(ijk_c(0), ijk_c(1), ijk_c(2));
+      const int k = marching_cubes_reoder[i];
+      codes(c,k) = code;
+#ifndef NDEBUG
+      Eigen::RowVector3i ijk_c_check;
+      code2ijk(code,ijk_c_check(0),ijk_c_check(1),ijk_c_check(2));
+      //printf("ijk2code(%d,%d,%d) = %ld\n",
+      //       ijk_c(0),ijk_c(1),ijk_c(2),code);
+      //printf("  code2ijk(%ld) = %d,%d,%d\n",
+      //    code, ijk_c_check(0),ijk_c_check(1),ijk_c_check(2));
+      // Check that the code is correct
+      assert(ijk_c_check(0) == ijk_c(0) && ijk_c_check(1) == ijk_c(1) &&
+             ijk_c_check(2) == ijk_c(2) &&
+             "ijk2code and code2ijk are not consistent");
+#endif
+    }
+  }
+  // igl::unique is the bottleneck by far.
+  // unique_rows might actually be faster because it doesn't do roundtrips to
+  // std::vector
+  Eigen::VectorXi I;
+  {
+    Eigen::Matrix<std::int64_t,Eigen::Dynamic,1> _;
+    Eigen::Matrix<typename DerivedJ::Scalar,Eigen::Dynamic,1> Jvec;
+    // igl::unique has a lot of internal copies :-(
+    igl::unique(codes,_,I,Jvec);
+    // reshape Jvec into ijk.rows() by 8
+    J = Jvec.reshaped(ijk.rows(),8);
+  }
+
+  // Map of codes as vector
+  const auto codes_vec = codes.reshaped();
+
+  unique_ijk.resize(I.size(),3);
+  for(int c = 0;c<I.size();c++)
+  {
+    code2ijk(codes_vec(I(c)),unique_ijk(c,0),unique_ijk(c,1),unique_ijk(c,2));
+  }
+}
+
+template<
+  typename Derivedorigin,
+  typename Derivedijk,
+  typename Derivedunique_ijk,
+  typename DerivedJ,
+  typename Derivedunique_corners
+    >
+IGL_INLINE void igl::unique_sparse_voxel_corners(
+  const Eigen::MatrixBase<Derivedorigin> & origin,
+  const typename Derivedorigin::Scalar h0,
+  const int depth,
+  const Eigen::MatrixBase<Derivedijk> & ijk,
+  Eigen::PlainObjectBase<Derivedunique_ijk> & unique_ijk,
+  Eigen::PlainObjectBase<DerivedJ> & J,
+  Eigen::PlainObjectBase<Derivedunique_corners> & unique_corners)
+{
+  unique_sparse_voxel_corners(depth,ijk,unique_ijk,J);
+
+  using Scalar = typename Derivedunique_corners::Scalar;
+  const Scalar h = h0 / (1 << depth);
+  unique_corners.resize(unique_ijk.rows(),3);
+  for(int c = 0;c<unique_ijk.rows();c++)
+  {
+    unique_corners(c,0) = origin(0) + h * unique_ijk(c,1);
+    unique_corners(c,1) = origin(1) + h * unique_ijk(c,0);
+    unique_corners(c,2) = origin(2) + h * unique_ijk(c,2);
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::unique_sparse_voxel_corners<Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 8, 1, -1, 8>>(int, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3>> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3>>&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 8, 1, -1, 8>>&);
+template void igl::unique_sparse_voxel_corners<Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 8, 1, -1, 8>, Eigen::Matrix<double, -1, 3, 1, -1, 3>>(Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3>> const&, Eigen::Matrix<double, 1, 3, 1, 1, 3>::Scalar, int, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3>> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3>>&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 8, 1, -1, 8>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>>&);
+#endif
+

include/igl/unique_sparse_voxel_corners.h

@@ -0,0 +1,64 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2025 Alec Jacobson <[email protected]>
+//
+// 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/.
+#ifndef IGL_UNIQUE_SPARSE_VOXEL_CORNERS_H
+#define IGL_UNIQUE_SPARSE_VOXEL_CORNERS_H
+
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl 
+{
+  ///  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] 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
+  template<
+    typename Derivedijk,
+    typename Derivedunique_ijk,
+    typename DerivedJ
+      >
+  IGL_INLINE void unique_sparse_voxel_corners(
+    const int depth,
+    const Eigen::MatrixBase<Derivedijk> & ijk,
+    Eigen::PlainObjectBase<Derivedunique_ijk> & unique_ijk,
+    Eigen::PlainObjectBase<DerivedJ> & J);
+
+  /// \brief overload
+  ///
+  ///   @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[out] unique_corners #unique_ijk by 3 list of unique corner
+  ///     positions
+  template<
+    typename Derivedorigin,
+    typename Derivedijk,
+    typename Derivedunique_ijk,
+    typename DerivedJ,
+    typename Derivedunique_corners
+      >
+  IGL_INLINE void unique_sparse_voxel_corners(
+    const Eigen::MatrixBase<Derivedorigin> & origin,
+    const typename Derivedorigin::Scalar h0,
+    const int depth,
+    const Eigen::MatrixBase<Derivedijk> & ijk,
+    Eigen::PlainObjectBase<Derivedunique_ijk> & unique_ijk,
+    Eigen::PlainObjectBase<DerivedJ> & J,
+    Eigen::PlainObjectBase<Derivedunique_corners> & unique_corners);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#    include "unique_sparse_voxel_corners.cpp"
+#endif
+#endif

tutorial/1001_LipschitzOctree/main.cpp

@@ -0,0 +1,237 @@
+#include <igl/read_triangle_mesh.h>
+#include <igl/lipschitz_octree.h>
+#include <igl/opengl/glfw/Viewer.h>
+#include <igl/AABB.h>
+#include <igl/fast_winding_number.h>
+#include <igl/signed_distance.h>
+#include <igl/grid.h>
+#include <igl/colormap.h>
+#include <igl/matlab_format.h>
+#include <igl/marching_cubes.h>
+#include <igl/unique_sparse_voxel_corners.h>
+#include <igl/get_seconds.h>
+
+
+// Edges of a cube box with minimum corner at origin and side length h.
+void box_edges(
+  const Eigen::RowVector3d & origin,
+  const double h,
+  Eigen::MatrixXd & V,
+  Eigen::MatrixXi & E)
+{
+  V.resize(8,3);
+  V<< 0,0,0,  0,0,1,  0,1,0,  0,1,1,  1,0,0,  1,0,1,  1,1,0,  1,1,1;
+  V *= h;
+  V.rowwise() += origin;
+  E.resize(12,2);
+  E<< 0,1,  0,2,  0,4,  1,3,  1,5,  2,3,  2,6,  3,7,  4,5,  4,6,  5,7,  6,7;
+}
+
+// Edges of cubes placed at all corners of octree with origin and cell side length h.
+void all_box_edges(
+  const Eigen::RowVector3d & origin,
+  const double h,
+  const Eigen::Matrix<int,Eigen::Dynamic,3,Eigen::RowMajor> & ijk,
+  Eigen::MatrixXd & V,
+  Eigen::MatrixXi & E)
+{
+  V.resize(ijk.rows()*8,3);
+  E.resize(ijk.rows()*12,2);
+  for(int c = 0;c<ijk.rows();c++)
+  {
+    Eigen::MatrixXd Vc;
+    Eigen::MatrixXi Ec;
+    box_edges(
+      origin + h * Eigen::RowVector3d(ijk(c,1),ijk(c,0),ijk(c,2)), h, Vc,Ec);
+    V.middleRows(c*8,8) = Vc;
+    E.middleRows(c*12,12) = Ec.array() + c*8;
+  }
+}
+
+int main(int argc, char * argv[])
+{
+  IGL_TICTOC_LAMBDA;
+  // Read mesh
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  if(!igl::read_triangle_mesh
+     (argc>1?argv[1]: TUTORIAL_SHARED_PATH "/decimated-knight.off",V,F)) {
+    std::cout << "Failed to load mesh." << std::endl;
+  }
+
+  // Prepare unsigned and signed distance function handles
+  igl::AABB<Eigen::MatrixXd,3> aabb;
+  aabb.init(V,F);
+  const double bbd = (V.colwise().maxCoeff() - V.colwise().minCoeff()).norm();
+  double offset_factor = 0.05;
+  igl::FastWindingNumberBVH fwn_bvh;
+  igl::fast_winding_number(V.cast<float>().eval(), F, 2, fwn_bvh);
+  const std::function<double(const Eigen::RowVector3d &)>
+    sdf = [&]( const Eigen::RowVector3d & p) -> double
+  {
+    double w = fast_winding_number(fwn_bvh, 2, p.template cast<float>().eval());         
+    double s = 1.-2.*std::abs(w);  
+    int i;
+    Eigen::RowVector3d c;
+    return s*std::sqrt(aabb.squared_distance(V,F,p,i,c)) - offset_factor*bbd;
+  };
+  // This can't be auto when using static library or compiler gets confused.
+  const std::function<double(const Eigen::RowVector3d &)>
+    udf = [&]( const Eigen::RowVector3d & p) -> double
+  {
+    return std::abs(sdf(p));
+  };
+
+  // Centered bounding cube (root of octree) with padding.
+  double h0 = (V.colwise().maxCoeff() - V.colwise().minCoeff()).maxCoeff();
+  int max_depth = 5;
+  // Pad bounding box by max_depth
+  //{
+  //  double leaf_h = h0 / (1 << max_depth);
+  //  h0 += leaf_h*2;
+  //}
+  // Pad to root
+  h0 *= 2;
+  Eigen::RowVector3d origin = 0.5*(V.colwise().maxCoeff() + V.colwise().minCoeff()).array() - 0.5*h0;
+
+
+  igl::opengl::glfw::Viewer viewer;
+  viewer.data().set_mesh(V,F);
+  const int main_index = viewer.selected_data_index;
+  Eigen::MatrixXd rV;
+  Eigen::MatrixXi rE;
+  box_edges(origin,h0,rV,rE);
+  viewer.data().set_edges(rV,rE,Eigen::RowVector3d(1,1,0));
+  viewer.data().set_face_based(true);
+  viewer.data().show_lines = false;
+  viewer.append_mesh();
+  const int aux_index = viewer.selected_data_index;
+
+  const auto update = [&]()
+  {
+    if(max_depth<=8)
+    {
+      // For comparison, how long would it take to compute on the dense grid
+      tictoc();
+      Eigen::Matrix<double,Eigen::Dynamic,3,Eigen::RowMajor> GV;
+      Eigen::RowVector3i res(1<<max_depth,1<<max_depth,1<<max_depth);
+      igl::grid(res,GV);
+      GV *= h0;
+      GV.rowwise() += origin;
+      Eigen::VectorXd GS(GV.rows());
+      igl::parallel_for( GV.rows(), [&](const int u)
+      {
+        // evaluate the function at the corner
+        GS(u) = sdf(GV.row(u));
+      },1000);
+      printf("                    n³ grid: %0.7f seconds\n",tictoc());
+    }else
+    {
+      printf("                    n³ grid: [omitted]\n");
+    }
+
+    tictoc();
+    // Identify the octree cells that could contain the surface
+    Eigen::Matrix<int,Eigen::Dynamic,3,Eigen::RowMajor> ijk;
+    igl::lipschitz_octree(origin,h0,max_depth,udf,ijk);
+    printf("           lipschitz_octree: %0.7f seconds\n",tictoc());
+
+    // Gather the corners of those leaf cells
+    const double h = h0 / (1 << max_depth);
+    Eigen::Matrix<int,Eigen::Dynamic,8,Eigen::RowMajor> J;
+    Eigen::Matrix<int,Eigen::Dynamic,3,Eigen::RowMajor> unique_ijk;
+    Eigen::Matrix<double,Eigen::Dynamic,3,Eigen::RowMajor> unique_corner_positions;
+    igl::unique_sparse_voxel_corners(origin,h0,max_depth,ijk,unique_ijk,J,unique_corner_positions);
+    printf("unique_sparse_voxel_corners: %0.7f seconds\n",tictoc());
+    /// Evaluate the signed distance function at the corners
+    Eigen::VectorXd S(unique_corner_positions.rows());
+    //for(int u = 0;u<unique_corner_positions.rows();u++)
+    igl::parallel_for(
+      unique_corner_positions.rows(),
+      [&](const int u)
+      {
+        // evaluate the function at the corner
+        S(u) = sdf(unique_corner_positions.row(u));
+      },1000);
+      printf("                        sdf: %0.7f seconds\n",tictoc());
+    // Run marching cubes on the sparse set of leaf cells
+    Eigen::Matrix<double,Eigen::Dynamic,3> mV;
+    Eigen::Matrix<int,Eigen::Dynamic,3> mF;
+    igl::marching_cubes( S,unique_corner_positions,J, 0, mV,mF);
+    printf("             marching_cubes: %0.7f seconds\n",tictoc());
+
+    // Visualize
+    viewer.data(aux_index).clear();
+    viewer.data(aux_index).set_mesh(mV,mF);
+    viewer.data(aux_index).set_face_based(true);
+    viewer.data().show_lines = max_depth <= 7;
+    viewer.data(aux_index).set_colors(Eigen::RowVector3d(0.941176,0.305882,0.282353));
+    if(max_depth <= 7)
+    {
+      // If tree is shallow enough, visualize the edges of the octree cells
+      Eigen::MatrixXd gV;
+      Eigen::MatrixXi gE;
+      all_box_edges(origin,h,ijk,gV,gE);
+      viewer.data(aux_index).set_edges(gV,gE,Eigen::RowVector3d(1,1,1));
+      Eigen::MatrixXd point_colors;
+      igl::colormap(
+          igl::ColorMapType::COLOR_MAP_TYPE_TURBO,
+          S,-S.array().abs().maxCoeff(),S.array().abs().maxCoeff(),
+          point_colors);
+      viewer.data(aux_index).add_points(unique_corner_positions,point_colors);
+      viewer.data(aux_index).point_size = 8;
+    }
+    printf("                     |ijk| : %ld\n",ijk.rows());
+  };
+
+  viewer.callback_key_pressed =
+  [&](igl::opengl::glfw::Viewer & viewer, unsigned int key, int mod)->bool
+  {
+    switch(key) {
+      default:
+        return false;
+      case '=':
+      case '+':
+      case '-':
+      case '_':
+        max_depth = std::max(0,max_depth + ((key == '-'||key=='_') ? -1 : 1));
+        printf("max_depth = %d\n",max_depth);
+        update();
+        break;
+      case '}':
+      case '{':
+        offset_factor =
+          std::min(std::max(-1.0,offset_factor+(key=='{'?-0.01:0.01)),1.0);
+        printf("offset_factor = %f\n",offset_factor);
+        update();
+        break;
+    }
+    return true;
+  };
+  update();
+  // Trivial batched function
+  const std::function<Eigen::VectorXd(const Eigen::Matrix<double,Eigen::Dynamic,3,Eigen::RowMajor> &)>
+    udf_batch = [&](const Eigen::Matrix<double,Eigen::Dynamic,3,Eigen::RowMajor> & P) -> Eigen::VectorXd
+  {
+    Eigen::VectorXd U(P.rows());
+    for(int i = 0;i<P.rows();i++)
+    {
+      U(i) = udf(P.row(i));
+    }
+    return U;
+  };
+  {
+    // The batched version should produce the same size output
+    Eigen::Matrix<int,Eigen::Dynamic,3,Eigen::RowMajor> ijk_batch;
+    igl::lipschitz_octree<true>(origin,h0,max_depth,udf_batch,ijk_batch);
+    printf("               |ijk_batch| : %ld\n",ijk_batch.rows());
+  }
+
+
+  std::cout << R"(Usage:
+  -,+  Decrease,Increase the depth of the octree and recompute
+  {,}  Decrease,Increase the offset factor for the signed distance function
+)";
+  viewer.launch();
+  return 0;
+}

tutorial/CMakeLists.txt

@@ -8,6 +8,7 @@ option(LIBIGL_TUTORIALS_CHAPTER6 "Compile libigl tutorials chapter 6" ON)
 option(LIBIGL_TUTORIALS_CHAPTER7 "Compile libigl tutorials chapter 7" ON)
 option(LIBIGL_TUTORIALS_CHAPTER8 "Compile libigl tutorials chapter 8" ON)
 option(LIBIGL_TUTORIALS_CHAPTER9 "Compile libigl tutorials chapter 9" ON)
+option(LIBIGL_TUTORIALS_CHAPTER10 "Compile libigl tutorials chapter 10" ON)
 
 # Download tutorial data
 include(libigl_tutorial_data)
@@ -131,3 +132,7 @@ if(LIBIGL_TUTORIALS_CHAPTER9)
     igl_add_tutorial(909_BatchMarchingCubes igl::glfw)
     igl_add_tutorial(910_OrientedBoundingBox igl::glfw igl_copyleft::cgal)
 endif()
+
+if(LIBIGL_TUTORIALS_CHAPTER10)
+  igl_add_tutorial(1001_LipschitzOctree igl::glfw)
+endif()