fix bugs and document better

include/igl/sparse_voxel_grid.cpp

@@ -59,16 +59,38 @@ IGL_INLINE void igl::sparse_voxel_grid(const Eigen::MatrixBase<DerivedP0>& p0,
   {
     Eigen::RowVector3i pi = queue.back();
     queue.pop_back();
+    if(visited.count(pi)){ continue; }
 
     VertexRowVector ctr = p0 + eps*pi.cast<ScalarV>(); // R^3 center of this cube
 
     // X, Y, Z basis vectors, and array of neighbor offsets used to construct cubes
     const Eigen::RowVector3i bx(1, 0, 0), by(0, 1, 0), bz(0, 0, -1);
-    const std::array<Eigen::RowVector3i, 6> neighbors = {
-      bx, -bx, by, -by, bz, -bz
-    };
+    const std::array<Eigen::RowVector3i, 26> neighbors = {
+      bx, -bx, by, -by, bz, -bz,
+      by-bz, -by+bz, // 1-2 4-7
+      bx+by, -bx-by, // 0-1 7-6
+      by+bz, -by-bz, // 0-3 6-5
+      by-bx, -by+bx, // 2-3 5-4
+      bx-bz, -bx+bz, // 1-5 3-7
+      bx+bz, -bx-bz, // 0-4 2-6
+      -bx+by+bz, bx-by-bz, // 3 5
+      bx+by+bz, -bx-by-bz, // 0 6
+      bx+by-bz, -bx-by+bz, //1 7
+      -bx+by-bz, bx-by+bz, // 2 4,
+  };
 
     // Compute the position of the cube corners and the scalar values at those corners
+    //
+    // Cube corners are ordered y-x-z, so their xyz offsets are:
+    //
+    // +++
+    // ++-
+    // -+-
+    // -++
+    // +-+
+    // +--
+    // ---
+    // --+
     std::array<VertexRowVector, 8> cubeCorners = {
       ctr+half_eps*(bx+by+bz).cast<ScalarV>(), ctr+half_eps*(bx+by-bz).cast<ScalarV>(), ctr+half_eps*(-bx+by-bz).cast<ScalarV>(), ctr+half_eps*(-bx+by+bz).cast<ScalarV>(),
       ctr+half_eps*(bx-by+bz).cast<ScalarV>(), ctr+half_eps*(bx-by-bz).cast<ScalarV>(), ctr+half_eps*(-bx-by-bz).cast<ScalarV>(), ctr+half_eps*(-bx-by+bz).cast<ScalarV>()
@@ -92,23 +114,51 @@ IGL_INLINE void igl::sparse_voxel_grid(const Eigen::MatrixBase<DerivedP0>& p0,
     // Add the cube vertices and indices to the output arrays if they are not there already
     IndexRowVector cube;
     uint8_t vertexAlreadyAdded = 0; // This is a bimask. If a bit is 1, it has been visited already by the BFS
-    constexpr std::array<uint8_t, 6> zv = {
+    constexpr std::array<uint8_t, 26> zv = {
       (1 << 0) | (1 << 1) | (1 << 4) | (1 << 5),
       (1 << 2) | (1 << 3) | (1 << 6) | (1 << 7),
       (1 << 0) | (1 << 1) | (1 << 2) | (1 << 3),
       (1 << 4) | (1 << 5) | (1 << 6) | (1 << 7),
       (1 << 0) | (1 << 3) | (1 << 4) | (1 << 7),
-      (1 << 1) | (1 << 2) | (1 << 5) | (1 << 6), };
-    constexpr std::array<std::array<int, 4>, 6> zvv {{
-        {{0, 1, 4, 5}}, {{3, 2, 7, 6}}, {{0, 1, 2, 3}},
-        {{4, 5, 6, 7}}, {{0, 3, 4, 7}}, {{1, 2, 5, 6}} }};
-
-    for (int n = 0; n < 6; n++) { // For each neighbor, check the hash table to see if its been added before
+      (1 << 1) | (1 << 2) | (1 << 5) | (1 << 6),
+      (1 << 1) | (1 << 2),
+      (1 << 4) | (1 << 7),
+      (1 << 0) | (1 << 1),
+      (1 << 6) | (1 << 7),
+      (1 << 0) | (1 << 3),
+      (1 << 5) | (1 << 6),
+      (1 << 2) | (1 << 3),
+      (1 << 4) | (1 << 5),
+      (1 << 1) | (1 << 5),
+      (1 << 3) | (1 << 7),
+      (1 << 0) | (1 << 4),
+      (1 << 2) | (1 << 6),
+      (1 << 3), (1 << 5), // diagonals
+      (1 << 0), (1 << 6),
+      (1 << 1), (1 << 7),
+      (1 << 2), (1 << 4),
+    };
+    constexpr std::array<std::array<int, 4>, 26> zvv {{
+      {{0, 1, 4, 5}}, {{3, 2, 7, 6}}, {{0, 1, 2, 3}},
+      {{4, 5, 6, 7}}, {{0, 3, 4, 7}}, {{1, 2, 5, 6}},
+      {{-1,-1,1,2}}, {{-1,-1,4,7}}, {{-1,-1,0,1}},{{-1,-1,7,6}},
+      {{-1,-1,0,3}}, {{-1,-1,5,6}}, {{-1,-1,2,3}}, {{-1,-1,5,4}},
+      {{-1,-1,1,5}}, {{-1,-1,3,7}}, {{-1,-1,0,4}}, {{-1,-1,2,6}},
+      {{-1,-1,-1,3}}, {{-1,-1,-1,5}}, {{-1,-1,-1,0}}, {{-1,-1,-1,6}},
+      {{-1,-1,-1,1}}, {{-1,-1,-1,7}}, {{-1,-1,-1,2}}, {{-1,-1,-1,4}} }};
+
+    for (int n = 0; n < 26; n++) { // For each neighbor, check the hash table to see if its been added before
       Eigen::RowVector3i nkey = pi + neighbors[n];
       auto nbr = visited.find(nkey);
       if (nbr != visited.end()) { // We've already visited this neighbor, use references to its vertices instead of duplicating them
         vertexAlreadyAdded |= zv[n];
-        for (int i = 0; i < 4; i++) { cube[zvv[n][i]] = CI_vector[nbr->second][zvv[n % 2 == 0 ? n + 1 : n - 1][i]]; }
+        for (int i = 0; i < 4; i++) 
+        {
+          if (zvv[n][i]!=-1)
+          {
+            cube[zvv[n][i]] = CI_vector[nbr->second][zvv[n % 2 == 0 ? n + 1 : n - 1][i]];
+          }
+        }
       } else {
         queue.push_back(nkey); // Otherwise, we have not visited the neighbor, put it in the BFS queue
       }
@@ -143,6 +193,9 @@ IGL_INLINE void igl::sparse_voxel_grid(const Eigen::MatrixBase<DerivedP0>& p0,
 
 
 #ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::sparse_voxel_grid<Eigen::Matrix<double, 1, 3, 1, 1, 3>, std::function<double (Eigen::Matrix<double, 1, 3, 1, 1, 3> const&)>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 8, 0, -1, 8> >(Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, std::function<double (Eigen::Matrix<double, 1, 3, 1, 1, 3> const&)> const&, double, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 8, 0, -1, 8> >&);
 template void igl::sparse_voxel_grid<class Eigen::Matrix<double, -1, -1, 0, -1, -1>, class std::function<double(class Eigen::Matrix<double, -1, -1, 0, -1, -1> const &)>, class Eigen::Matrix<double, -1, 1, 0, -1, 1>, 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 std::function<double(class Eigen::Matrix<double, -1, -1, 0, -1, -1> const &)> const &, double, int, class Eigen::PlainObjectBase<class Eigen::Matrix<double, -1, 1, 0, -1, 1> > &, class Eigen::PlainObjectBase<class Eigen::Matrix<double, -1, -1, 0, -1, -1> > &, class Eigen::PlainObjectBase<class Eigen::Matrix<int, -1, -1, 0, -1, -1> > &);
 template void igl::sparse_voxel_grid<Eigen::Matrix<double, 1, 3, 1, 1, 3>, std::function<double (Eigen::Matrix<double, 1, 3, 1, 1, 3> const&)>, 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::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, std::function<double (Eigen::Matrix<double, 1, 3, 1, 1, 3> const&)> const&, double, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::sparse_voxel_grid<Eigen::Matrix<double, 1, 3, 1, 1, 3>, std::function<double (Eigen::Matrix<double, 1, 3, 1, 1, 3> const&)>, 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::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, std::function<double (Eigen::Matrix<double, 1, 3, 1, 1, 3> const&)> const&, double, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);

include/igl/sparse_voxel_grid.h

@@ -9,11 +9,10 @@
 #define IGL_SPARSE_VOXEL_GRID_H
 
 #include "igl_inline.h"
-
 #include <Eigen/Core>
 
-namespace igl {
-
+namespace igl 
+{
   // sparse_voxel_grid( p0, scalarFunc, eps, CV, CS, CI )
   //
   // Given a point, p0, on an isosurface, construct a shell of epsilon sized cubes surrounding the surface.
@@ -27,22 +26,30 @@ namespace igl {
   //  eps  The edge length of the cubes surrounding the surface
   //  expected_number_of_cubes  This pre-allocates internal data structures to speed things up
   // Output:
-  //   CS  #cube-vertices by 1 list of scalar values at the cube vertices
-  //   CV  #cube-vertices by 3 list of cube vertex positions
-  //   CI  #number of cubes by 8 list of indexes into CS and CV. Each row represents a cube
+  //   CS  #CV by 1 list of scalar values at the cube vertices
+  //   CV  #CV by 3 list of cube vertex positions
+  //   CI  #CI by 8 list of cube indices into rows of CS and CV. Each row
+  //     represents 8 corners of cube in y-x-z binary counting order.
   //
-  template <typename DerivedP0, typename Func, typename DerivedS, typename DerivedV, typename DerivedI>
-  IGL_INLINE void sparse_voxel_grid(const Eigen::MatrixBase<DerivedP0>& p0,
-                                    const Func& scalarFunc,
-                                    const double eps,
-                                    const int expected_number_of_cubes,
-                                    Eigen::PlainObjectBase<DerivedS>& CS,
-                                    Eigen::PlainObjectBase<DerivedV>& CV,
-                                    Eigen::PlainObjectBase<DerivedI>& CI);
+  template <
+    typename DerivedP0, 
+    typename Func, 
+    typename DerivedS, 
+    typename DerivedV, 
+    typename DerivedI>
+  IGL_INLINE void sparse_voxel_grid(
+    const Eigen::MatrixBase<DerivedP0>& p0,
+    const Func& scalarFunc,
+    const double eps,
+    const int expected_number_of_cubes,
+    Eigen::PlainObjectBase<DerivedS>& CS,
+    Eigen::PlainObjectBase<DerivedV>& CV,
+    Eigen::PlainObjectBase<DerivedI>& CI);
 
 }
+
 #ifndef IGL_STATIC_LIBRARY
 #    include "sparse_voxel_grid.cpp"
 #endif
 
-#endif // IGL_SPARSE_VOXEL_GRID_H
+#endif