Minimal AABB tree + SDFs + Variable Radius Offsets (#2490)

* working eytzinger aabb

* working and reasonably efficient variable radius offset

* working example

* comments

* better key commands

* bad includes, defines

* missing sign function; test

include/igl/SphereMeshWedge.cpp

@@ -0,0 +1,231 @@
+// 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 "SphereMeshWedge.h"
+#include "round_cone_signed_distance.h"
+#include "sign.h"
+#include <Eigen/QR>
+#include <Eigen/Geometry>
+
+template <typename Scalar>
+IGL_INLINE igl::SphereMeshWedge<Scalar>::SphereMeshWedge(
+  const RowVector3S & V0,
+  const RowVector3S & V1,
+  const RowVector3S & V2,
+  const Scalar r0,
+  const Scalar r1,
+  const Scalar r2)
+{
+  // Internal copy
+  V.row(0) = V0;
+  V.row(1) = V1;
+  V.row(2) = V2;
+  r(0) = r0;
+  r(1) = r1;
+  r(2) = r2;
+
+  flavor = FULL;
+  // By default use full
+
+  EV.row(0) = V.row(2) - V.row(1);
+  EV.row(1) = V.row(0) - V.row(2);
+  EV.row(2) = V.row(1) - V.row(0);
+  l = EV.rowwise().norm();
+  l2 = l.array().square();
+  rr << r(1) - r(2), r(2) - r(0), r(0) - r(1);
+  a2 = l2.array() - rr.array().square();
+  il2 = 1.0/l2.array();
+
+  /////////////////////////////////////////////
+  /// BIG_VERTEX ?
+  /////////////////////////////////////////////
+  {
+    r.maxCoeff(&max_i);
+    int j = (max_i+1)%3;
+    int k = (max_i+2)%3;
+    if((l(k) + r(j) < r(max_i)) && (l(j) + r(k) < r(max_i)))
+    {
+      flavor = BIG_VERTEX;
+    }
+  }
+
+  /////////////////////////////////////////////
+  /// BIG_EDGE ?
+  /////////////////////////////////////////////
+  if(flavor == FULL)
+  {
+    // Case where one edge's roundCone containes the others
+    for(int e = 0;e<3;e++)
+    {
+      const int i = (e+1)%3;
+      const int j = (e+2)%3;
+      const int k = (e+3)%3;
+      const Scalar s = 
+        igl::round_cone_signed_distance(V.row(k),V.row(i),V.row(j),r(i),r(j));
+      if(-s > r(k))
+      {
+        flavor = BIG_EDGE;
+        max_i = i;
+        break;
+      }
+    }
+  }
+
+  if(flavor == FULL && !compute_planes())
+  {
+    flavor = NO_TRIANGLE;
+  }
+}
+
+template <typename Scalar>
+IGL_INLINE Scalar igl::SphereMeshWedge<Scalar>::operator()(const RowVector3S & p) const
+{
+  if(flavor == BIG_VERTEX)
+  {
+    // Case 0: Vertex i
+    return (p - V.row(max_i)).norm() - r(max_i);
+  }
+
+  if(flavor == BIG_EDGE)
+  {
+    const int i = max_i;
+    const int j = (i+1)%3;
+    // Case 1: Edge e
+    return this->round_cone_signed_distance(p,i,j);
+  }
+
+  Scalar s = std::numeric_limits<Scalar>::infinity();
+  if(flavor == FULL)
+  {
+    // This is possibly the bottleneck and could be turned into precomputed
+    // plane equations.
+
+    // signed distance to triangle plane (this is immediately recomputed later in
+    // sdSkewedExtrudedTriangle...)
+    const auto plane_sdf = [](
+        const RowVector3S & p,
+        const Eigen::RowVector4d & plane)
+    {
+      return plane.head<3>().dot(p) + plane(3);
+    };
+    Scalar d0 = plane_sdf(p, planes.row(0));
+    Scalar planes_s = -std::abs(d0);
+    // Reflect if necessary so that q is always on negative side of plane
+    RowVector3S q = p - (d0 - planes_s) * planes.row(0).template head<3>();
+    // Other planes (for negative side slab, by symmetry)
+    for(int i = 1;i<planes.rows();i++)
+    {
+      planes_s = std::max(planes_s,plane_sdf(q, planes.row(i)));
+    }
+    // This produces correct interior distance
+    if(planes_s <= 0)
+    {
+      s = std::min(s,planes_s);
+    }else
+    {
+
+      const auto & nor = planes.row(1).template head<3>();
+      const RowVector3S q0 = q - T.row(0);
+      const RowVector3S q1 = q - T.row(1);
+      const RowVector3S q2 = q - T.row(2);
+
+      if(!(sign(C.row(0).dot(q0)) + 
+           sign(C.row(1).dot(q1)) + 
+           sign(C.row(2).dot(q2))<2.0))
+      {
+        s = std::min(s,planes_s);
+      }
+    }
+
+    //s = std::min(s,sdSkewedExtrudedTriangle(q,V,T));
+    //s = std::min(s,sdSkewedExtrudedTriangle(p,B,V));
+  }
+  assert(flavor == FULL || flavor == NO_TRIANGLE);
+
+  for(int e = 0;e<3;e++)
+  {
+    const int i = (e+1)%3;
+    const int j = (e+2)%3;
+    s = std::min(s,this->round_cone_signed_distance(p,i,j));
+  }
+  return s;
+}
+
+template <typename Scalar>
+IGL_INLINE bool igl::SphereMeshWedge<Scalar>::compute_planes()
+{
+  // Non-degenerate case
+  const RowVector3S & a = V.row(0);
+  const RowVector3S & b = V.row(1);
+  const RowVector3S & c = V.row(2);
+  const Scalar & ra = r(0);
+  const Scalar & rb = r(1);
+  const Scalar & rc = r(2);
+  Eigen::Matrix<Scalar,2,3,Eigen::RowMajor> A;
+  A<<
+    b-a,
+    c-a;
+  const Eigen::Vector2d d(rb-ra,rc-ra);
+  const RowVector3S N = (A.row(0).cross(A.row(1))).normalized();
+  //const Eigen::CompleteOrthogonalDecomposition<decltype(A)> cod(A);
+  const RowVector3S n0 = A.completeOrthogonalDecomposition().solve(d);
+  const Scalar qA = N.squaredNorm();
+  // qB is zeros by construction. We could delete all terms involving qB
+  // It's not even clear if keeping them would lead to more accurate results.
+  const Scalar qB = 2 * N.dot(n0);
+  const Scalar qC = n0.squaredNorm() - 1;
+  const Scalar qD = qB*qB - 4*qA*qC;
+
+  if(qD<0) { return false; }
+
+  Scalar t_sol_1 = (-qB + std::sqrt(qD)) / (2*qA);
+  RowVector3S n1 = -(t_sol_1 * N + n0);
+  T  = V + r * n1;
+
+  const auto plane_equation = [](
+      const RowVector3S & a,
+      const RowVector3S & b,
+      const RowVector3S & c)->Eigen::RowVector4d
+  {
+    RowVector3S n = (b-a).cross(c-a).normalized();
+    n.normalize();
+    Scalar d = -n.dot(a);
+    return Eigen::RowVector4d(n(0),n(1),n(2),d);
+  };
+
+  planes.row(0) = plane_equation(V.row(0),V.row(1),V.row(2));
+  planes.row(1) = plane_equation(T.row(2),T.row(1),T.row(0));
+  planes.row(2) = plane_equation(V.row(1),V.row(0),T.row(0));
+  planes.row(3) = plane_equation(V.row(2),V.row(1),T.row(1));
+  planes.row(4) = plane_equation(V.row(0),V.row(2),T.row(2));
+
+  // Determine if the closest point is on the face.
+  const RowVector3S v10 = T.row(1) - T.row(0); 
+  const RowVector3S v21 = T.row(2) - T.row(1); 
+  const RowVector3S v02 = T.row(0) - T.row(2); 
+  const auto & nor = planes.row(1).template head<3>();
+  const RowVector3S c10 = v10.cross(nor);
+  const RowVector3S c21 = v21.cross(nor);
+  const RowVector3S c02 = v02.cross(nor);
+  C<<c10,c21,c02;
+  return true;
+}
+
+template <typename Scalar>
+IGL_INLINE Scalar igl::SphereMeshWedge<Scalar>::round_cone_signed_distance(const RowVector3S & p, const int i, const int j) const
+{
+  const int e = (j+1)%3;
+  return igl::round_cone_signed_distance(
+    p, V.row(i), r(i), r(j), EV.row(e), l2(e), rr(e), a2(e), il2(e));
+}
+
+
+#ifdef IGL_STATIC_LIBRARY
+/// Explicit template instantiation
+template class igl::SphereMeshWedge<double>;
+#endif
+

include/igl/SphereMeshWedge.h

@@ -0,0 +1,87 @@
+// 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_SPHERE_MESH_WEDGE_H
+#define IGL_SPHERE_MESH_WEDGE_H
+
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  /// A class to compute the signed distance to a "Sphere-Mesh Wedge" as seen in
+  /// variable radius offset surfaces or Sphere-Meshes. Each wedge is defined
+  /// by three vertices and three radii, one at each vertex. The wedge is
+  /// the union of all spheres at points on the triangle with radius linearly
+  /// interpolated. See, e.g., "Sphere-Meshes for Real-Time Hand Modeling and
+  /// Tracking" or "A Multilinear Model for Bidirectional Craniofacial
+  /// Reconstruction" or "Sphere-Meshes: Shape Approximation using Spherical
+  /// Quadric Error Metrics" or "Variable-Radius Offset Surface Approximation on
+  /// the GPU".
+  ///
+  template <typename Scalar>
+  class SphereMeshWedge
+  {
+    public: 
+      using RowVector3S = Eigen::Matrix<Scalar, 1, 3>;
+      // Fields
+      enum 
+      {
+        BIG_VERTEX = 0,
+        BIG_EDGE = 1,
+        NO_TRIANGLE = 2,
+        FULL = 3
+      } flavor;
+      Eigen::Matrix<Scalar,3,3,Eigen::RowMajor> V;
+      Eigen::Matrix<Scalar,3,1> r;
+      Eigen::Matrix<Scalar,3,3,Eigen::RowMajor> EV;
+      Eigen::Matrix<Scalar,3,1> l,l2,rr,a2,il2;
+      int max_i;
+      Eigen::Matrix<Scalar,5,4,Eigen::RowMajor> planes;
+      Eigen::Matrix<Scalar,3,3,Eigen::RowMajor> T;
+      Eigen::Matrix<Scalar,3,3,Eigen::RowMajor> C;
+
+      SphereMeshWedge(){}
+      /// Constructor that takes three vertices and three radii
+      ///
+      /// @param V0  first vertex position
+      /// @param V1  second vertex position
+      /// @param V2  third vertex position
+      /// @param r0  radius at first vertex
+      /// @param r1  radius at second vertex
+      /// @param r2  radius at third vertex
+      IGL_INLINE SphereMeshWedge(
+        const RowVector3S & V0,
+        const RowVector3S & V1,
+        const RowVector3S & V2,
+        const Scalar r0,
+        const Scalar r1,
+        const Scalar r2);
+      /// @param[in] p  3-vector query point
+      /// @return signed distance to the wedge at point p
+      IGL_INLINE Scalar operator()(const RowVector3S & p) const;
+    private:
+      /// Precompute planes used for determining bounded signed to the skewed
+      /// triangular slab portion.
+      ///
+      /// @return true if planes are well defined (false implies this slab has
+      /// no contribution).
+      IGL_INLINE bool compute_planes();
+      /// Compute the signed distance to the wedge at a point p for the edge
+      /// (i,j) 
+      ///
+      /// @param[in] p  3-vector query point
+      /// @param[in] i  index of first vertex (0,1,2)
+      /// @param[in] j  index of second vertex (0,1,2)
+      IGL_INLINE Scalar round_cone_signed_distance(const RowVector3S & p, const int i, const int j) const;
+  };
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#include "SphereMeshWedge.cpp"
+#endif
+#endif 

include/igl/adjacency_matrix.cpp

@@ -111,6 +111,8 @@ IGL_INLINE void igl::adjacency_matrix(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::adjacency_matrix<Eigen::Matrix<int, -1, 3, 1, -1, 3>, int>(Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3>> const&, Eigen::SparseMatrix<int, 0, int>&);
 template void igl::adjacency_matrix<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, int>(Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::SparseMatrix<int, 0, int>& );
 // generated by autoexplicit.sh
 template void igl::adjacency_matrix<Eigen::Matrix<int, -1, -1, 0, -1, -1>, bool>(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<bool, 0, int>&);

include/igl/box_faces.h

@@ -1,3 +1,10 @@
+// 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_BOX_FACES_H
 #define IGL_BOX_FACES_H
 #include "igl_inline.h"

include/igl/box_simplices.cpp

@@ -0,0 +1,43 @@
+// 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 "box_simplices.h"
+
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedB
+>
+IGL_INLINE void igl::box_simplices(
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
+  Eigen::PlainObjectBase<DerivedB> & B1,
+  Eigen::PlainObjectBase<DerivedB> & B2)
+{
+  B1.setConstant(F.rows(),V.cols(),std::numeric_limits<double>::infinity());
+  B2.setConstant(F.rows(),V.cols(),-std::numeric_limits<double>::infinity());
+  for(int f = 0; f < F.rows(); f++) 
+  {
+    for(int c = 0; c < F.cols(); c++) 
+    {
+      for(int d = 0; d < V.cols(); d++) 
+      {
+        B1(f,d) = std::min(B1(f,d),V(F(f,c),d));
+        B2(f,d) = std::max(B2(f,d),V(F(f,c),d));
+      }
+    }
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::box_simplices<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 2, 1, -1, 2>, Eigen::Matrix<double, -1, 3, 1, -1, 3>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 2, 1, -1, 2>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>>&);
+// generated by autoexplicit.sh
+template void igl::box_simplices<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<double, -1, 3, 1, -1, 3>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>>&);
+template void igl::box_simplices<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3>>&);
+#endif 

include/igl/box_simplices.h

@@ -0,0 +1,37 @@
+// 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_BOX_SIMPLICES_H
+#define IGL_BOX_SIMPLICES_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  /// 
+  /// 
+  /// @param[in] V  #V by dim list of vertex positions
+  /// @param[in] F #F by ss list of simplex indices into V
+  /// @param[out] B1  #B by dim list of minimum corners of the Eytzinger AABBs
+  /// @param[out] B2  #B by dim list of maximum corners of the Eytzinger AABBs
+  ///
+  template <
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedB
+  >
+  IGL_INLINE void box_simplices(
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    Eigen::PlainObjectBase<DerivedB> & B1,
+    Eigen::PlainObjectBase<DerivedB> & B2);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "box_simplices.cpp"
+#endif
+#endif 

include/igl/edges.cpp

@@ -64,6 +64,8 @@ IGL_INLINE void igl::edges(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::edges<Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 2, 1, -1, 2>>(Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3>> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 1, -1, 2>>&);
+// generated by autoexplicit.sh
 template void igl::edges<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::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::edges<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 2, 0, -1, 2> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&);
 template void igl::edges<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);

include/igl/eytzinger_aabb.cpp

@@ -0,0 +1,90 @@
+#include "eytzinger_aabb.h"
+#include "PlainMatrix.h"
+#include "median.h"
+#include <vector>
+#include <algorithm>
+
+template <
+  typename DerivedPB,
+  typename DerivedB,
+  typename Derivedleaf
+>
+IGL_INLINE void igl::eytzinger_aabb(
+  const Eigen::MatrixBase<DerivedPB> & PB1,
+  const Eigen::MatrixBase<DerivedPB> & PB2,
+  Eigen::PlainObjectBase<DerivedB> & B1,
+  Eigen::PlainObjectBase<DerivedB> & B2,
+  Eigen::PlainObjectBase<Derivedleaf> & leaf)
+{
+  using Scalar = typename DerivedPB::Scalar;
+  using VectorXS = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
+  // Number of primitives
+  const int m = PB1.rows();
+  assert(PB1.rows() == PB2.rows());
+  if(m==0)
+  {
+    B1.resize(0,PB1.cols());
+    B2.resize(0,PB2.cols());
+    leaf.resize(0);
+    return;
+  }
+  const PlainMatrix<DerivedPB> PBC = (PB1 + PB2) / 2;
+  const int complete_m = 1 << ((int)std::ceil(std::log2(m))+1);
+  B1.resize(complete_m, PB1.cols());
+  B2.resize(complete_m, PB2.cols());
+  leaf.resize(complete_m);
+  leaf.setConstant(-2);
+  std::vector<int> I(m); for(int i = 0; i < m; i++) { I[i] = i; }
+
+  std::function<void(const int, const std::vector<int> &)> recursive_helper;
+  recursive_helper = [&](const int i, const std::vector<int> & I)
+  {
+    if(I.size() == 0) { return; }
+
+    B1.row(i) = PB1(I, Eigen::all).colwise().minCoeff();
+    B2.row(i) = PB2(I, Eigen::all).colwise().maxCoeff();
+
+    if(I.size() == 1) { leaf(i) = I[0]; return; }
+
+    leaf(i) = -1;
+    int dir;
+    (B2.row(i) - B1.row(i)).maxCoeff(&dir);
+
+    Scalar split_value;
+    VectorXS Z = PBC(I, dir);
+    igl::median(Z,split_value); 
+    std::vector<int> left; left.reserve(I.size()/2+1);
+    std::vector<int> right; right.reserve(I.size()/2+1);
+    for(int j = 0; j < (int)I.size(); j++)
+    {
+      if(Z(j) < split_value)
+      {
+        left.push_back(I[j]);
+      }else if(Z(j) > split_value)
+      {
+        right.push_back(I[j]);
+      }
+    }
+    for(int j = 0; j < (int)I.size(); j++)
+    {
+      if(Z(j) == split_value)
+      {
+        (left.size()<right.size() ? left : right).push_back(I[j]);
+      }
+    }
+    assert(std::abs(int(left.size()-right.size())) <= 1);
+
+    const int left_i = 2*i + 1;
+    const int right_i = 2*i + 2;
+    recursive_helper(left_i, left);
+    recursive_helper(right_i, right);
+  };
+  recursive_helper(0,I);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::eytzinger_aabb<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::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>>&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>>&);
+template void igl::eytzinger_aabb<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3>>&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>>&);
+#endif

include/igl/eytzinger_aabb.h

@@ -0,0 +1,40 @@
+// 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_EYTZINGER_AABB_H
+#define IGL_EYTZINGER_AABB_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  /// Compute the Eytzinger AABB for a given mesh
+  /// 
+  /// @param[in] PB1  #P by dim list of minimum corners of the AABBs
+  /// @param[in] PB2  #P by dim list of maximum corners of the AABBs
+  /// @param[out] B1  #B by dim list of minimum corners of the Eytzinger AABBs
+  /// @param[out] B2  #B by dim list of maximum corners of the Eytzinger AABBs
+  /// @param[out] leaf #B list of leaf indices, -1 indicates internal node, -2
+  /// indicates empty node
+  ///
+  template <
+    typename DerivedPB,
+    typename DerivedB,
+    typename Derivedleaf
+  >
+  IGL_INLINE void eytzinger_aabb(
+    const Eigen::MatrixBase<DerivedPB> & PB1,
+    const Eigen::MatrixBase<DerivedPB> & PB2,
+    Eigen::PlainObjectBase<DerivedB> & B1,
+    Eigen::PlainObjectBase<DerivedB> & B2,
+    Eigen::PlainObjectBase<Derivedleaf> & leaf);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "eytzinger_aabb.cpp"
+#endif
+#endif 

include/igl/eytzinger_aabb_sdf.cpp

@@ -0,0 +1,100 @@
+#include "eytzinger_aabb_sdf.h"
+#include "parallel_for.h"
+
+template <
+  typename Derivedp,
+  typename Func,
+  typename DerivedB,
+  typename Derivedleaf
+>
+IGL_INLINE void igl::eytzinger_aabb_sdf(
+  const Eigen::MatrixBase<Derivedp> & p,
+  const Func & primitive,
+  const Eigen::MatrixBase<DerivedB> & B1,
+  const Eigen::MatrixBase<DerivedB> & B2,
+  const Eigen::MatrixBase<Derivedleaf> & leaf,
+  typename Derivedp::Scalar & f)
+{
+  using Scalar = typename Derivedp::Scalar;
+  using RowVectorS = Eigen::Matrix<Scalar, 1, Derivedp::ColsAtCompileTime>;
+  // https://iquilezles.org/articles/distfunctions/
+  const auto box_sdf = [&p](const RowVectorS& b1, const RowVectorS& b2) -> Scalar
+  {
+    // half diagonal
+    const RowVectorS b = 0.5*(b2 - b1);
+    // relative position
+    const RowVectorS r = p - 0.5*(b1 + b2);
+    const RowVectorS q = r.cwiseAbs() - b;
+    const Scalar t1 = q.cwiseMax(0.0).squaredNorm();
+    const Scalar t2 = std::min(q.maxCoeff(),0.0);
+    return std::sqrt(t1 + t2);
+  };
+
+  // List of current active nodes
+  std::vector<std::pair<int,Scalar>> active;
+  active.reserve(50);
+  active.emplace_back(0,box_sdf(B1.row(0), B2.row(0)));
+  f = std::numeric_limits<Scalar>::infinity();
+  while(!active.empty())
+  {
+    const auto [i, box_f] = active.back();
+    active.pop_back();
+    if(box_f >= f){ continue; }
+    if(leaf(i) >= 0)
+    {
+      f = std::min(f,primitive(leaf(i)));
+    }else
+    {
+      const int left_i = 2*i + 1;
+      const int right_i = 2*i + 2;
+      const Scalar left_f = box_sdf(B1.row(left_i), B2.row(left_i));
+      const Scalar right_f = box_sdf(B1.row(right_i), B2.row(right_i));
+      if(left_f < right_f)
+      {
+        active.emplace_back(right_i, right_f);
+        active.emplace_back(left_i, left_f);
+      }else
+      {
+        active.emplace_back(left_i, left_f);
+        active.emplace_back(right_i, right_f);
+      }
+    }
+  }
+}
+
+template <
+  typename DerivedP,
+  typename Func,
+  typename DerivedB,
+  typename Derivedleaf,
+  typename DerivedS
+>
+IGL_INLINE void igl::eytzinger_aabb_sdf(
+  const Eigen::MatrixBase<DerivedP> & P,
+  const Func & primitive,
+  const Eigen::MatrixBase<DerivedB> & B1,
+  const Eigen::MatrixBase<DerivedB> & B2,
+  const Eigen::MatrixBase<Derivedleaf> & leaf,
+  Eigen::PlainObjectBase<DerivedS> & S)
+{
+  using Scalar = typename DerivedS::Scalar;
+  S.resize(P.rows());
+  //for(int i = 0; i < P.rows(); i++)
+  igl::parallel_for(P.rows(), [&](const int i)
+  {
+    const Eigen::Matrix<Scalar,1,DerivedP::ColsAtCompileTime> p = P.row(i);
+    const std::function<Scalar(const int)> primitive_i = [&](const int j)
+    {
+      return primitive(p,j);
+    };
+    igl::eytzinger_aabb_sdf(p,primitive_i,B1,B2,leaf,S(i));
+  },1000);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::eytzinger_aabb_sdf<Eigen::Matrix<double, -1, 3, 1, -1, 3>, std::function<double (Eigen::Matrix<double, 1, 3, 1, 1, 3> const&, int)>, Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>> const&, std::function<double (Eigen::Matrix<double, 1, 3, 1, 1, 3> const&, int)> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1>>&);
+template void igl::eytzinger_aabb_sdf<Eigen::Matrix<double, 1, 3, 1, 1, 3>, std::function<double (int)>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>>(Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3>> const&, std::function<double (int)> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>> const&, Eigen::Matrix<double, 1, 3, 1, 1, 3>::Scalar&);
+template void igl::eytzinger_aabb_sdf<Eigen::Matrix<double, -1, 3, 0, -1, 3>, std::function<double (Eigen::Matrix<double, 1, 3, 1, 1, 3> const&, int)>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3>> const&, std::function<double (Eigen::Matrix<double, 1, 3, 1, 1, 3> const&, int)> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1>>&);
+#endif

include/igl/eytzinger_aabb_sdf.h

@@ -0,0 +1,66 @@
+// 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_EYTZINGER_AABB_SDF_H
+#define IGL_EYTZINGER_AABB_SDF_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  /// 
+  /// @param[in] B1  #B by dim list of minimum corners of the Eytzinger AABBs
+  /// @param[in] B2  #B by dim list of maximum corners of the Eytzinger AABBs
+  /// @param[in] leaf #B list of leaf indices, -1 indicates internal node, -2
+  /// indicates empty node
+  /// @param[in] primitive  function handle that takes as input a primitive id
+  /// and returns the SDF to that primitive `primitive(i)`
+  /// @param[out] f  SDF value at query point x
+  ///
+  template <
+    typename Derivedp,
+    typename Func,
+    typename DerivedB,
+    typename Derivedleaf
+  >
+  IGL_INLINE void eytzinger_aabb_sdf(
+    const Eigen::MatrixBase<Derivedp> & p,
+    const Func & primitive,
+    const Eigen::MatrixBase<DerivedB> & B1,
+    const Eigen::MatrixBase<DerivedB> & B2,
+    const Eigen::MatrixBase<Derivedleaf> & leaf,
+    typename Derivedp::Scalar & f);
+
+  /// @param[in] P  #P by dim list of query points
+  /// @param[in] B1  #B by dim list of minimum corners of the Eytzinger AABBs
+  /// @param[in] B2  #B by dim list of maximum corners of the Eytzinger AABBs
+  /// @param[in] leaf #B list of leaf indices, -1 indicates internal node, -2
+  /// indicates empty node
+  /// @param[in] primitive  function handle that takes as input the query point
+  ///   (row of P) and a primitive id and returns the SDF to that primitive
+  ///   `primitive(p,i)`
+  /// @param[out] S  #P list of SDF values at query points P
+  template <
+    typename DerivedP,
+    typename Func,
+    typename DerivedB,
+    typename Derivedleaf,
+    typename DerivedS
+  >
+  IGL_INLINE void eytzinger_aabb_sdf(
+    const Eigen::MatrixBase<DerivedP> & P,
+    const Func & primitive,
+    const Eigen::MatrixBase<DerivedB> & B1,
+    const Eigen::MatrixBase<DerivedB> & B2,
+    const Eigen::MatrixBase<Derivedleaf> & leaf,
+    Eigen::PlainObjectBase<DerivedS> & S);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "eytzinger_aabb_sdf.cpp"
+#endif
+#endif 

include/igl/grid.cpp

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

include/igl/march_cube.cpp

@@ -134,6 +134,9 @@ 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, unsigned int, 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<unsigned int, 8, 1, 0, 8, 1> const&, double const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3>>&, unsigned int&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -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, 3, 0, -1, 3>>, double, unsigned int, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3>> 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, 0, -1, 3>>&, unsigned int&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -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>>>&);
+// 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> > >&);

include/igl/marching_cubes.cpp

@@ -149,6 +149,9 @@ IGL_INLINE void igl::marching_cubes(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::marching_cubes<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 3, 1, -1, 3>, 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&, unsigned int, unsigned int, unsigned int, 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>>&);
+template void igl::marching_cubes<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, 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, 0, -1, 3>> const&, unsigned int, unsigned int, unsigned int, 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>>&);
 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> >&);

include/igl/matlab_format.cpp

@@ -8,6 +8,8 @@
 #include "matlab_format.h"
 #include "STR.h"
 #include "find.h"
+#include "list_to_matrix.h"
+
 
 template <typename DerivedM>
 IGL_INLINE const Eigen::WithFormat< DerivedM > igl::matlab_format(
@@ -42,6 +44,17 @@ IGL_INLINE std::string igl::matlab_format_index(
   return STR(igl::matlab_format((M.array()+1).eval(),name));
 }
 
+template <typename T>
+IGL_INLINE std::string igl::matlab_format_index(
+  const std::vector<T> & Mvec,
+  const std::string name)
+{
+  Eigen::Matrix<T,Eigen::Dynamic,1> M;
+  igl::list_to_matrix(Mvec,M);
+  // can't return WithFormat since that uses a pointer to matrix
+  return STR(igl::matlab_format((M.array()+1).eval(),name));
+}
+
 template <typename DerivedS>
 IGL_INLINE const std::string
 igl::matlab_format(
@@ -118,6 +131,7 @@ IGL_INLINE const std::string igl::matlab_format(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+template std::basic_string<char, std::char_traits<char>, std::allocator<char>> igl::matlab_format_index<int>(std::vector<int, std::allocator<int>> const&, std::basic_string<char, std::char_traits<char>, std::allocator<char>>);
 // generated by autoexplicit.sh
 template std::basic_string<char, std::char_traits<char>, std::allocator<char> > const igl::matlab_format<bool>(Eigen::SparseMatrix<bool, 0, int> const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
 // generated by autoexplicit.sh

include/igl/matlab_format.h

@@ -13,6 +13,7 @@
 #include <Eigen/Core>
 #include <Eigen/Sparse>
 #include <string>
+#include <vector>
 
 namespace igl
 {
@@ -51,6 +52,10 @@ namespace igl
   IGL_INLINE std::string matlab_format_index(
     const Eigen::MatrixBase<DerivedM> & M,
     const std::string name = "");
+  template <typename T>
+  IGL_INLINE std::string matlab_format_index(
+    const std::vector<T> & Mvec,
+    const std::string name = "");
   /// Same but for sparse matrices. Print IJV format into an auxiliary variable
   /// and then print a call to sparse which will construct the sparse matrix
   ///

include/igl/round_cone_signed_distance.cpp

@@ -0,0 +1,70 @@
+// 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 "round_cone_signed_distance.h"
+#include "sign.h"
+
+template <
+  typename Derivedp,
+  typename Deriveda, 
+  typename Derivedb>
+IGL_INLINE typename Derivedp::Scalar igl::round_cone_signed_distance(
+  const Eigen::MatrixBase<Derivedp> & p, 
+  const Eigen::MatrixBase<Deriveda> & a, 
+  const Eigen::MatrixBase<Derivedb> & b, 
+  const typename Derivedp::Scalar & r1, 
+  const typename Derivedp::Scalar & r2)
+{
+  // https://iquilezles.org/articles/distfunctions/
+  using Scalar = typename Derivedp::Scalar;
+  using RowVector3S = Eigen::Matrix<Scalar, 1, 3>;
+  // sampling independent computations (only depend on shape)
+  const RowVector3S ba = b - a;
+  const Scalar l2 = ba.dot(ba);
+  const Scalar rr = r1 - r2;
+  const Scalar a2 = l2 - rr*rr;
+  const Scalar il2 = 1.0/l2;
+  return round_cone_signed_distance(p,a,r1,r2,ba,l2,rr,a2,il2);
+}
+
+template <
+  typename Derivedp,
+  typename Deriveda, 
+  typename Derivedba>
+IGL_INLINE typename Derivedp::Scalar igl::round_cone_signed_distance(
+  const Eigen::MatrixBase<Derivedp> & p, 
+  const Eigen::MatrixBase<Deriveda> & a, 
+  const typename Derivedp::Scalar & r1, 
+  const typename Derivedp::Scalar & r2,
+  const Eigen::MatrixBase<Derivedba> & ba,
+  const typename Derivedp::Scalar & l2,
+  const typename Derivedp::Scalar & rr,
+  const typename Derivedp::Scalar & a2,
+  const typename Derivedp::Scalar & il2)
+{
+  using Scalar = typename Derivedp::Scalar;
+  using RowVector3S = Eigen::Matrix<Scalar, 1, 3>;
+  // sampling dependant computations
+  const RowVector3S pa = p - a;
+  const Scalar y = pa.dot(ba);
+  const Scalar z = y - l2;
+  const Scalar x2 = ( pa*l2 - ba*y ).squaredNorm();
+  const Scalar y2 = y*y*l2;
+  const Scalar z2 = z*z*l2;
+
+  // single square root!
+  const Scalar k = sign(rr)*rr*rr*x2;
+  if( sign(z)*a2*z2>k ) return  sqrt(x2 + z2)        *il2 - r2;
+  if( sign(y)*a2*y2<k ) return  sqrt(x2 + y2)        *il2 - r1;
+                        return (sqrt(x2*a2*il2)+y*rr)*il2 - r1;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template Eigen::Block<Eigen::Matrix<double, 3, 3, 1, 3, 3>, 1, 3, true>::Scalar igl::round_cone_signed_distance<Eigen::Block<Eigen::Matrix<double, 3, 3, 1, 3, 3>, 1, 3, true>, Eigen::Block<Eigen::Matrix<double, 3, 3, 1, 3, 3>, 1, 3, true>, Eigen::Block<Eigen::Matrix<double, 3, 3, 1, 3, 3>, 1, 3, true>>(Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, 3, 3, 1, 3, 3>, 1, 3, true>> const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, 3, 3, 1, 3, 3>, 1, 3, true>> const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, 3, 3, 1, 3, 3>, 1, 3, true>> const&, Eigen::Block<Eigen::Matrix<double, 3, 3, 1, 3, 3>, 1, 3, true>::Scalar const&, Eigen::Block<Eigen::Matrix<double, 3, 3, 1, 3, 3>, 1, 3, true>::Scalar const&);
+template Eigen::Matrix<double, 1, 3, 1, 1, 3>::Scalar igl::round_cone_signed_distance<Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Block<Eigen::Matrix<double, 3, 3, 1, 3, 3> const, 1, 3, true>, Eigen::Block<Eigen::Matrix<double, 3, 3, 1, 3, 3> const, 1, 3, true>>(Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3>> const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, 3, 3, 1, 3, 3> const, 1, 3, true>> const&, Eigen::Matrix<double, 1, 3, 1, 1, 3>::Scalar const&, Eigen::Matrix<double, 1, 3, 1, 1, 3>::Scalar const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, 3, 3, 1, 3, 3> const, 1, 3, true>> const&, Eigen::Matrix<double, 1, 3, 1, 1, 3>::Scalar const&, Eigen::Matrix<double, 1, 3, 1, 1, 3>::Scalar const&, Eigen::Matrix<double, 1, 3, 1, 1, 3>::Scalar const&, Eigen::Matrix<double, 1, 3, 1, 1, 3>::Scalar const&);
+#endif

include/igl/round_cone_signed_distance.h

@@ -0,0 +1,66 @@
+// 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_ROUND_CONE_SIGNED_DISTANCE_H
+#define IGL_ROUND_CONE_SIGNED_DISTANCE_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  /// Compute signed distance to a round cone (sphere-mesh capsule). 
+  ///
+  /// @param[in] p  3-vector query point
+  /// @param[in] a  3-vector position of the first vertex
+  /// @param[in] b  3-vector position of the second vertex
+  /// @param[in] r1  radius at vertex a
+  /// @param[in] r2  radius at vertex b
+  /// @return signed distance to the round cone
+  template <
+    typename Derivedp,
+    typename Deriveda, 
+    typename Derivedb>
+  IGL_INLINE typename Derivedp::Scalar round_cone_signed_distance(
+    const Eigen::MatrixBase<Derivedp> & p, 
+    const Eigen::MatrixBase<Deriveda> & a, 
+    const Eigen::MatrixBase<Derivedb> & b, 
+    const typename Derivedp::Scalar & r1, 
+    const typename Derivedp::Scalar & r2);
+
+  /// Compute signed distance to a round cone (sphere-mesh capsule) with
+  /// pre-computed values that are independent of the query point.
+  ///
+  /// @param[in] p  3-vector query point
+  /// @param[in] a  3-vector position of the first vertex
+  /// @param[in] r1  radius at vertex a
+  /// @param[in] r2  radius at vertex b
+  /// @param[in] ba  3-vector direction from a to b
+  /// @param[in] l2  squared length of ba
+  /// @param[in] rr  r1 - r2
+  /// @param[in] a2  l2 - rr*rr
+  /// @param[in] il2  1/l2
+  /// @return signed distance to the round cone
+  template <
+    typename Derivedp,
+    typename Deriveda, 
+    typename Derivedba>
+  IGL_INLINE typename Derivedp::Scalar round_cone_signed_distance(
+    const Eigen::MatrixBase<Derivedp> & p, 
+    const Eigen::MatrixBase<Deriveda> & a, 
+    const typename Derivedp::Scalar & r1, 
+    const typename Derivedp::Scalar & r2,
+    const Eigen::MatrixBase<Derivedba> & ba,
+    const typename Derivedp::Scalar & l2,
+    const typename Derivedp::Scalar & rr,
+    const typename Derivedp::Scalar & a2,
+    const typename Derivedp::Scalar & il2);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "round_cone_signed_distance.cpp"
+#endif
+#endif 

include/igl/sign.h

@@ -0,0 +1,20 @@
+// 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_SIGN_H
+#define IGL_SIGN_H
+namespace igl
+{
+  /// Returns the sign of a number. https://en.wikipedia.org/wiki/Sign_function
+  ///
+  /// @param[in] x The number to compute the sign of.
+  /// @return 1 if x > 0, -1 if x < 0, and 0 if x == 0.
+  ///
+  /// https://stackoverflow.com/questions/1903954/is-there-a-standard-sign-function-signum-sgn-in-c-c
+  template <typename T> inline T sign(const T & x){return (x>T(0)) - (x<T(0));}
+}
+#endif

include/igl/variable_radius_offset.cpp

@@ -0,0 +1,195 @@
+// 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 "variable_radius_offset.h"
+#include "SphereMeshWedge.h"
+#include "marching_cubes.h"
+#include "unique_sparse_voxel_corners.h"
+#include "lipschitz_octree.h"
+#include "eytzinger_aabb.h"
+#include "eytzinger_aabb_sdf.h"
+#include "get_seconds.h"
+#include "parallel_for.h"
+#include <cstdio>
+
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedR,
+  typename Derivedorigin,
+  typename DerivedmV,
+  typename DerivedmF>
+void igl::variable_radius_offset(
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    const Eigen::MatrixBase<DerivedR> & R,
+    const Eigen::MatrixBase<Derivedorigin> & origin,
+    const typename Derivedorigin::Scalar h0,
+    const int max_depth,
+    Eigen::PlainObjectBase<DerivedmV> & mV,
+    Eigen::PlainObjectBase<DerivedmF> & mF)
+{
+  using Scalar = typename DerivedV::Scalar;
+  Eigen::Matrix<Scalar,Eigen::Dynamic,3,Eigen::RowMajor> PB1,PB2;
+  std::vector<igl::SphereMeshWedge<Scalar>> data;
+  std::function<Scalar(const Eigen::Matrix<Scalar,1,3> &,const int i)> primitive;
+  igl::variable_radius_offset(V,F,R,PB1,PB2,data,primitive);
+  igl::variable_radius_offset(
+    PB1,PB2,data,primitive,origin,h0,max_depth,mV,mF);
+}
+
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedR,
+  typename DerivedPB1,
+  typename DerivedPB2,
+  typename Scalar>
+void igl::variable_radius_offset(
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
+  const Eigen::MatrixBase<DerivedR> & R,
+  Eigen::PlainObjectBase<DerivedPB1> & PB1,
+  Eigen::PlainObjectBase<DerivedPB2> & PB2,
+  std::vector<igl::SphereMeshWedge<Scalar>> & data,
+  std::function<Scalar(const Eigen::Matrix<Scalar,1,3> &,const int i)> & 
+    primitive)
+{
+  // Bounds
+  PB1.setConstant(F.rows(),3,std::numeric_limits<Scalar>::infinity());
+  PB2.setConstant(F.rows(),3,-std::numeric_limits<Scalar>::infinity());
+  for(int f = 0;f<F.rows();f++)
+  {
+    for(int c = 0;c<F.cols();c++)
+    {
+      for(int d = 0;d<V.cols();d++)
+      {
+        PB1(f,d) = std::min(PB1(f,d),V(F(f,c),d)-R(F(f,c)));
+        PB2(f,d) = std::max(PB2(f,d),V(F(f,c),d)+R(F(f,c)));
+      }
+    }
+  }
+
+  // Precomputed data
+  data.resize(F.rows());
+  for(int f = 0;f<F.rows();f++)
+  {
+    data[f] = igl::SphereMeshWedge<Scalar>(
+      V.row(F(f,0)),V.row(F(f,1)),V.row(F(f,2)),
+      R(F(f,0)),R(F(f,1)),R(F(f,2)));
+  }
+
+  // I suppose primitive could just capture `data` by value. Then it owns the
+  // data and we don't need to worry about `data` getting destroyed before its
+  // called. I'm not sure how expensive that copy would be.
+  primitive = [&](const Eigen::RowVector3d & p, const int i)
+  {
+    return data[i](p);
+  };
+}
+
+
+template <
+  typename DerivedPB1,
+  typename DerivedPB2,
+  typename Scalar,
+  typename Derivedorigin,
+  typename DerivedmV,
+  typename DerivedmF>
+void igl::variable_radius_offset(
+  const Eigen::MatrixBase<DerivedPB1> & PB1,
+  const Eigen::MatrixBase<DerivedPB2> & PB2,
+  const std::vector<igl::SphereMeshWedge<Scalar>> & data,
+  const std::function<Scalar(const Eigen::Matrix<Scalar,1,3> &,const int i)> & 
+    primitive,
+  const Eigen::MatrixBase<Derivedorigin> & origin,
+  const Scalar h0,
+  const int max_depth,
+  Eigen::PlainObjectBase<DerivedmV> & mV,
+  Eigen::PlainObjectBase<DerivedmF> & mF)
+{
+  using RowVector3S = Eigen::Matrix<Scalar,1,3>;
+  IGL_TICTOC_LAMBDA;
+  tictoc();
+  Eigen::Matrix<Scalar,Eigen::Dynamic,3,Eigen::RowMajor> B1,B2;
+  Eigen::VectorXi leaf;
+  igl::eytzinger_aabb(PB1,PB2,B1,B2,leaf);
+  printf("%-20s: %g secs\n","eytzinger_aabb",tictoc());
+
+  const std::function<Scalar(const RowVector3S &)>
+    sdf = [&](const RowVector3S & p) -> Scalar
+  {
+    const std::function<Scalar(const int)> primitive_p = [&](const int j)
+    {
+      return primitive(p,j);
+    };
+    Scalar f;
+    igl::eytzinger_aabb_sdf(p,primitive_p,B1,B2,leaf,f);
+    return f;
+  };
+  const std::function<Scalar(const RowVector3S &)>
+    udf = [&](const RowVector3S & p) -> Scalar
+  {
+    return std::abs(sdf(p));
+    //// This made performance worse.
+    //const std::function<Scalar(const int)> primitive_p = [&](const int j)
+    //{
+    //  const Scalar d = udTriangle( V.row(F(j,0)), V.row(F(j,1)), V.row(F(j,2)), p);
+    //  const RowVector3S r(R(F(j,0)),R(F(j,1)),R(F(j,2)));
+    //  const Scalar min_r = r.minCoeff();
+    //  const Scalar max_r = r.maxCoeff();
+    //  if(d > max_r)
+    //  {
+    //    return d - max_r;
+    //  }else if(d < min_r)
+    //  {
+    //    return d - min_r;
+    //  }
+    //  return 0.0;
+    //};
+    //Scalar f;
+    //igl::eytzinger_aabb_sdf(p,primitive_p,B1,B2,leaf,f);
+    //return f;
+  };
+  Eigen::Matrix<int,Eigen::Dynamic,3,Eigen::RowMajor> ijk;
+  igl::lipschitz_octree( origin,h0,max_depth,udf,ijk);
+  printf("%-20s: %g secs\n","lipschitz_octree",tictoc());
+
+  {
+    tictoc();
+    // Gather the corners of those leaf cells
+    const Scalar 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<Scalar,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());
+    printf("%-20s: %g secs\n","unique_sparse_vo...",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());
+      printf("%-20s: %g secs\n","sdf",tictoc());
+    // Run marching cubes on the sparse set of leaf cells
+    igl::marching_cubes( S,unique_corner_positions,J, 0, mV,mF);
+    //printf("             marching_cubes: %0.7f seconds\n",tictoc());
+    printf("%-20s: %g secs\n","marching_cubes",tictoc());
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::variable_radius_offset<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<double, -1, 3, 1, -1, 3>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>>&, std::vector<igl::SphereMeshWedge<double>, std::allocator<igl::SphereMeshWedge<double>>>&, std::function<double (Eigen::Matrix<double, 1, 3, 1, 1, 3> const&, int)>&);
+template void igl::variable_radius_offset<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, 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::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1>> const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3>> const&, Eigen::Matrix<double, 1, 3, 1, 1, 3>::Scalar, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3>>&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3>>&);
+#endif

include/igl/variable_radius_offset.h

@@ -0,0 +1,106 @@
+// 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_VARIABLE_RADIUS_OFFSET_H
+#define IGL_VARIABLE_RADIUS_OFFSET_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+#include <functional>
+
+namespace igl
+{
+  template <typename Scalar>
+  class SphereMeshWedge;
+
+  /// Compute the variable radius offset (see, e.g., "Variable-Radius Offset
+  /// Curves and Surfaces" [Lin & Rokne 1997]) of a triangle mesh with a
+  /// pieceiwse linear radius function.
+  ///
+  /// @param[in] V  #V by 3 list of vertex positions
+  /// @param[in] F  #F by 3 list of triangle indices into V
+  /// @param[in] R  #V list of radii for each vertex
+  /// @param[in] origin  3-vector origin of octree used for mesh extraction
+  /// @param[in] h0  side length of octree root cell
+  /// @param[in] max_depth  maximum depth of octree used for mesh extraction
+  ///   (root is depth=0)
+  /// @param[out] mV #mV by 3 list of vertex positions of the offset mesh
+  /// @param[out] mF #mF by 3 list of triangle indices into mV
+  template <
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedR,
+    typename Derivedorigin,
+    typename DerivedmV,
+    typename DerivedmF>
+  void variable_radius_offset(
+      const Eigen::MatrixBase<DerivedV> & V,
+      const Eigen::MatrixBase<DerivedF> & F,
+      const Eigen::MatrixBase<DerivedR> & R,
+      const Eigen::MatrixBase<Derivedorigin> & origin,
+      const typename Derivedorigin::Scalar h0,
+      const int max_depth,
+      Eigen::PlainObjectBase<DerivedmV> & mV,
+      Eigen::PlainObjectBase<DerivedmF> & mF);
+
+  /// \brief Overload which prepares data for a eytzinger_aabb_sdf callable sdf
+  /// function.
+  ///
+  /// @param[in] V  #V by 3 list of vertex positions
+  /// @param[in] F  #F by 3 list of triangle indices into V
+  /// @param[in] R  #V list of radii for each vertex
+  /// @param[out] PB1 #F by 3 list of minimum corners of each offset
+  ///   triangle's axis aligned bounding box
+  /// @param[out] PB2 #F by 3 list of maximum corners of each offset
+  ///   triangle's axis aligned bounding box
+  /// @param[out] #F list of precomputed sphere-mesh wedge data
+  /// @param[out] primitive  function handle that takes as input a point and
+  ///   primitive index as input and outputs the corresponding signed distance.
+  ///   This function assumes that `data` will remain alive.
+  template <
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedR,
+    typename DerivedPB1,
+    typename DerivedPB2,
+    typename Scalar>
+  void variable_radius_offset(
+      const Eigen::MatrixBase<DerivedV> & V,
+      const Eigen::MatrixBase<DerivedF> & F,
+      const Eigen::MatrixBase<DerivedR> & R,
+      Eigen::PlainObjectBase<DerivedPB1> & PB1,
+      Eigen::PlainObjectBase<DerivedPB2> & PB2,
+      std::vector<igl::SphereMeshWedge<Scalar>> & data,
+      std::function<Scalar(const Eigen::Matrix<Scalar,1,3> &,const int i)> & 
+        primitive);
+
+  ///// \brief Overload using precomputed data.
+  template <
+    typename DerivedPB1,
+    typename DerivedPB2,
+    typename Scalar,
+    typename Derivedorigin,
+    typename DerivedmV,
+    typename DerivedmF>
+  void variable_radius_offset(
+    const Eigen::MatrixBase<DerivedPB1> & PB1,
+    const Eigen::MatrixBase<DerivedPB2> & PB2,
+    const std::vector<igl::SphereMeshWedge<Scalar>> & data,
+    const std::function<Scalar(const Eigen::Matrix<Scalar,1,3> &,const int i)> & 
+      primitive,
+    const Eigen::MatrixBase<Derivedorigin> & origin,
+    const Scalar h0,
+    const int max_depth,
+    Eigen::PlainObjectBase<DerivedmV> & mV,
+    Eigen::PlainObjectBase<DerivedmF> & mF);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "variable_radius_offset.cpp"
+#endif
+#endif 
+

tests/include/igl/sign.cpp

@@ -0,0 +1,28 @@
+#include <test_common.h>
+#include <igl/sign.h>
+
+
+template <typename T>
+void test()
+{
+  REQUIRE (igl::sign(T(  0)) == 0);
+  REQUIRE (igl::sign(T( -0)) == 0);
+  REQUIRE (igl::sign(T(  1)) == 1);
+  REQUIRE (igl::sign(T( -1)) == -1);
+  if constexpr (std::is_floating_point<T>::value)
+  {
+    REQUIRE (igl::sign(T(  0.0)) == 0);
+    REQUIRE (igl::sign(T( -0.0)) == 0);
+    REQUIRE (igl::sign(T(0.5)) == 1);
+    REQUIRE (igl::sign(T(-0.5)) == -1);
+    REQUIRE (igl::sign(T(0.00000000000001)) == 1);
+    REQUIRE (igl::sign(T(-0.00000000000001)) == -1);
+  }
+}
+
+TEST_CASE("sign: cases", "[igl]" )
+{
+  test<float>();
+  test<double>();
+  test<int>();
+}

tutorial/1002_EytzingerAABB/main.cpp

@@ -0,0 +1,313 @@
+#include <igl/read_triangle_mesh.h>
+#include <igl/box_simplices.h>
+#include <igl/eytzinger_aabb.h>
+#include <igl/lipschitz_octree.h>
+#include <igl/unique_sparse_voxel_corners.h>
+#include <igl/eytzinger_aabb_sdf.h>
+#include <igl/unique.h>
+#include <igl/readDMAT.h>
+#include <igl/grid.h>
+#include <igl/edges.h>
+#include <igl/marching_cubes.h>
+#include <igl/parallel_for.h>
+#include <igl/opengl/glfw/Viewer.h>
+#include <igl/matlab_format.h>
+#include <igl/variable_radius_offset.h>
+#include <igl/get_seconds.h>
+#include <igl/SphereMeshWedge.h>
+
+#include <limits>
+
+#include <Eigen/Geometry>
+
+// What kind of offset is being computed?
+enum class OffsetType
+{
+  EDGE_OFFSET = 0,
+  TRIANGLE_OFFSET = 1,
+  VARIABLE_RADIUS_OFFSET = 2
+} offset = OffsetType::VARIABLE_RADIUS_OFFSET;
+
+///////////////////////////////////////////////////////////////////////
+/// Some helper functions for triangle and line segment distance
+///////////////////////////////////////////////////////////////////////
+template <typename T> inline T sign(const T & x){ return (x>T(0)) - (x<T(0)); }
+
+inline double dot2(const Eigen::RowVector3d & v)
+{
+  return v.dot(v);
+}
+
+inline double clamp(const double & x, const double & a, const double & b)
+{
+  return std::max(a, std::min(b, x));
+}
+
+// https://iquilezles.org/articles/triangledistance/
+template <typename T1, typename T2, typename T3, typename Tp>
+typename Tp::Scalar udTriangle(
+    const T1 & v1, const T2 & v2, const T3 & v3, const Tp & p)
+{
+  using vec3 = Eigen::Matrix<typename Tp::Scalar, 3, 1>;
+  // prepare data    
+  vec3 v21 = v2 - v1; vec3 p1 = p - v1;
+  vec3 v32 = v3 - v2; vec3 p2 = p - v2;
+  vec3 v13 = v1 - v3; vec3 p3 = p - v3;
+  vec3 nor =  v21.cross( v13 );
+
+  return sqrt( // inside/outside test    
+    (sign(v21.cross(nor).dot(p1)) + 
+     sign(v32.cross(nor).dot(p2)) + 
+     sign(v13.cross(nor).dot(p3))<2.0) 
+     ?
+     // 3 edges    
+     std::min( std::min( 
+     dot2(v21*clamp(v21.dot(p1)/dot2(v21),0.0,1.0)-p1), 
+     dot2(v32*clamp(v32.dot(p2)/dot2(v32),0.0,1.0)-p2) ), 
+     dot2(v13*clamp(v13.dot(p3)/dot2(v13),0.0,1.0)-p3) )
+     :
+     // 1 face    
+     nor.dot(p1)*nor.dot(p1)/dot2(nor) );
+}
+    
+template <typename Tp, typename Ta, typename Tb>
+typename Tp::Scalar udLineSegment(const Tp & p, const Ta & a, const Tb & b)
+{
+  using vec3 = Eigen::Matrix<typename Tp::Scalar, 3, 1>;
+  vec3 pa = p - a;
+  vec3 ba = b - a;
+  const auto h = clamp( pa.dot(ba)/ba.dot(ba), 0.0, 1.0 );
+  return ( pa - ba*h ).norm();
+}
+
+int main(int argc, char * argv[])
+{
+  IGL_TICTOC_LAMBDA;
+  // Octree depth
+  int max_depth = 7;
+
+  // Read mesh
+  Eigen::Matrix<double,Eigen::Dynamic,3,Eigen::RowMajor> V;
+  Eigen::Matrix<int,Eigen::Dynamic,3,Eigen::RowMajor> 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;
+  }
+  // Normalize
+  V.rowwise() -= 0.5* (V.colwise().maxCoeff() + V.colwise().minCoeff());
+  const double scale = 0.8/V.array().abs().maxCoeff();
+  V *= scale;
+
+  // Per-vertex radius for VARIABLE_RADIUS_OFFSET
+  // if second arg exists read R from .dmat
+  Eigen::VectorXd R;
+  if(argc>2)
+  {
+    igl::readDMAT( argv[2],R);
+    R *= scale;
+  }else
+  {
+    R = V.col(1);
+    R.array() -= R.minCoeff();
+    R /= R.maxCoeff();
+    R.array() += 1.0;
+    R = (R.array().pow(4)).eval();
+    R *= 0.01;
+  }
+
+  ////////////////////////////////////////////////////////////////////////
+  /// Box up simplices (and, for VARIABLE_RADIUS_OFFSET, precompute data)
+  ////////////////////////////////////////////////////////////////////////
+  tictoc();
+  Eigen::Matrix<double,Eigen::Dynamic,3,Eigen::RowMajor> PB1,PB2;
+  Eigen::Matrix<int,Eigen::Dynamic,2,Eigen::RowMajor> E;
+  igl::edges(F,E);
+
+  double isovalue = 0.02;
+  std::function<double(const Eigen::RowVector3d &,const int i)> primitive;
+  std::vector<igl::SphereMeshWedge<double>> data;
+  switch(offset)
+  {
+    case OffsetType::EDGE_OFFSET:
+    {
+      igl::box_simplices(V,E,PB1,PB2);
+      primitive = [&](const Eigen::RowVector3d & p, const int i)
+      {
+        return udLineSegment(p, V.row(E(i,0)), V.row(E(i,1)));
+      };
+      break;
+    }
+    case OffsetType::TRIANGLE_OFFSET:
+    {
+      igl::box_simplices(V,F,PB1,PB2);
+      primitive = [&](const Eigen::RowVector3d & p, const int i)
+      {
+        return udTriangle( V.row(F(i,0)), V.row(F(i,1)), V.row(F(i,2)), p);
+      };
+      break;
+    }
+    case OffsetType::VARIABLE_RADIUS_OFFSET:
+    {
+      isovalue = 0;
+      igl::variable_radius_offset(V,F,R,PB1,PB2,data,primitive);
+      break;
+    }
+  }
+  printf("%-20s: %g secs\n","PB",tictoc());
+
+  ////////////////////////////////////////////////////////////////////////
+  /// Put boxes into eytzinger AABB tree
+  ////////////////////////////////////////////////////////////////////////
+  tictoc();
+  Eigen::Matrix<double,Eigen::Dynamic,3,Eigen::RowMajor> B1,B2;
+  Eigen::VectorXi leaf;
+  igl::eytzinger_aabb(PB1,PB2,B1,B2,leaf);
+  printf("%-20s: %g secs\n","eytzinger_aabb",tictoc());
+  //{
+  //  Eigen::VectorXi U;
+  //  igl::unique(leaf,U);
+  //  printf("%d → %d\n",PB1.rows(),U.size()-2);
+  //}
+
+  igl::opengl::glfw::Viewer viewer;
+  Eigen::Matrix<double,Eigen::Dynamic,3> mV,oV;
+  Eigen::Matrix<int,Eigen::Dynamic,3> mF,oF;
+  viewer.data().set_mesh(V,F);
+  viewer.data().set_face_based(true);
+  viewer.data().show_lines = false;
+
+  if(offset == OffsetType::VARIABLE_RADIUS_OFFSET)
+  {
+    viewer.data().set_colors(R);
+  }
+  viewer.append_mesh();
+  const int m_index = viewer.selected_data_index;
+  viewer.append_mesh();
+  const int o_index = viewer.selected_data_index;
+
+  viewer.data(m_index).set_face_based(true);
+  viewer.data(m_index).show_lines = true;
+  viewer.data(o_index).set_face_based(true);
+  viewer.data(o_index).show_lines = false;
+
+
+  const auto update = [&]()
+  {
+    ////////////////////////////////////////////////////////////////////////
+    /// If the grid isn't that big then compute a dense M.C. mesh
+    /// Call batched eytzinger_aabb_sdf directly on grid nodes.
+    ////////////////////////////////////////////////////////////////////////
+    tictoc();
+    const int ns = (1 << max_depth) + 1;
+    if(ns <= 128+1)
+    {
+      Eigen::RowVector3i res(ns,ns,ns);
+      Eigen::Matrix<double,Eigen::Dynamic,3,Eigen::RowMajor> GV;
+      igl::grid(res,GV);
+      GV *= 2;
+      GV.array() -= 1;
+      Eigen::VectorXd S;
+      igl::eytzinger_aabb_sdf(GV,primitive,B1,B2,leaf,S);
+      printf("%-20s: %g secs\n","eytzinger_aabb_sdf",tictoc());
+
+      tictoc();
+      igl::marching_cubes(S,GV,res(0),res(1),res(2),isovalue,mV,mF);
+      printf("%-20s: %g secs\n","marching_cubes",tictoc());
+    }else
+    {
+      mV.resize(0,3); 
+      mF.resize(0,3);
+    }
+
+    ////////////////////////////////////////////////////////////////////////
+    /// Prepare an unsigned distance and signed distance function handle.
+    ////////////////////////////////////////////////////////////////////////
+    tictoc();
+    const std::function<double(const Eigen::RowVector3d &)>
+      sdf = [&](const Eigen::RowVector3d & p) -> double
+    {
+      const std::function<double(const int)> primitive_p = [&](const int j)
+      {
+        return primitive(p,j);
+      };
+      double f;
+      igl::eytzinger_aabb_sdf(p,primitive_p,B1,B2,leaf,f);
+      return f - isovalue;
+    };
+    const std::function<double(const Eigen::RowVector3d &)>
+      udf = [&](const Eigen::RowVector3d & p) -> double
+    {
+      return std::abs(sdf(p));
+    };
+    ////////////////////////////////////////////////////////////////////////
+    /// Use udf to build sparse voxel octree around the zero level set
+    ////////////////////////////////////////////////////////////////////////
+    Eigen::RowVector3d origin(-1,-1,-1);
+    const double h0 = 2;
+    Eigen::Matrix<int,Eigen::Dynamic,3,Eigen::RowMajor> ijk;
+    igl::lipschitz_octree( origin,h0,max_depth,udf,ijk);
+    printf("%-20s: %g secs\n","lipschitz_octree",tictoc());
+
+    ////////////////////////////////////////////////////////////////////////
+    /// Gather the corners of those leaf cells, compute sdf there and run 
+    /// (sparse) marching cubes
+    ////////////////////////////////////////////////////////////////////////
+    {
+      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());
+      printf("%-20s: %g secs\n","unique_sparse_vo...",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());
+        printf("%-20s: %g secs\n","sdf",tictoc());
+      // Run marching cubes on the sparse set of leaf cells
+      igl::marching_cubes( S,unique_corner_positions,J, 0, oV,oF);
+      //printf("             marching_cubes: %0.7f seconds\n",tictoc());
+      printf("%-20s: %g secs\n","marching_cubes",tictoc());
+    }
+
+    viewer.data(m_index).clear();
+    viewer.data(m_index).set_mesh(mV,mF);
+    viewer.data(m_index).set_colors(Eigen::RowVector3d(0.8,0.2,0.2));
+
+    viewer.data(o_index).clear();
+    viewer.data(o_index).set_mesh(oV,oF);
+    viewer.data(o_index).set_colors(Eigen::RowVector3d(0.2,0.8,0.2));
+  };
+
+  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;
+    }
+    return true;
+  };
+  update();
+
+  viewer.launch();
+  return 0;
+}

tutorial/910_OrientedBoundingBox/main.cpp

@@ -53,7 +53,10 @@ int main(int argc, char * argv[])
   Eigen::MatrixXi G;
   igl::copyleft::cgal::convex_hull(V,W,G);
   Eigen::Matrix3d BR;
+  igl::oriented_bounding_box(W, 50000, igl::ORIENTED_BOUNDING_BOX_MINIMIZE_SURFACE_AREA, BR);
+  std::cout<<igl::matlab_format(BR, "saBR")<<std::endl;
   igl::oriented_bounding_box(W, 50000, igl::ORIENTED_BOUNDING_BOX_MINIMIZE_VOLUME, BR);
+  std::cout<<igl::matlab_format(BR, "vBR")<<std::endl;
   double t_best = tictoc();
 
   // Cube mesh

tutorial/CMakeLists.txt

@@ -135,4 +135,5 @@ endif()
 
 if(LIBIGL_TUTORIALS_CHAPTER10)
   igl_add_tutorial(1001_LipschitzOctree igl::glfw)
+  igl_add_tutorial(1002_EytzingerAABB igl::glfw)
 endif()