Support Batched Marching Cubes (#2422)

* Working example; need to change name rather than overwrite 705

* separate tutorial for batch

* comment

include/igl/march_cube.cpp

@@ -37,91 +37,97 @@ IGL_INLINE void igl::march_cube(
   std::unordered_map<std::int64_t,int> & E2V)
 {
 
-// These consts get stored reasonably
+  // These consts get stored reasonably
 #include "marching_cubes_tables.h"
 
   // Seems this is also successfully inlined
+  //
+  // Returns whether the vertex is new
   const auto ij2vertex =
-    [&E2V,&V,&n,&GV]
-      (const Index & i, const Index & j, const Scalar & t)->Index
-  {
-    // Seems this is successfully inlined.
-    const auto ij2key = [](std::int32_t i,std::int32_t j)
+    [&E2V,&V,&n]
+    (const Index & i, const Index & j, Index & v)->bool
     {
-      if(i>j){ std::swap(i,j); }
-      std::int64_t ret = 0;
-      ret |= i;
-      ret |= static_cast<std::int64_t>(j) << 32;
-      return ret;
+      // Seems this is successfully inlined.
+      const auto ij2key = [](std::int32_t i,std::int32_t j)
+      {
+        if(i>j){ std::swap(i,j); }
+        std::int64_t ret = 0;
+        ret |= i;
+        ret |= static_cast<std::int64_t>(j) << 32;
+        return ret;
+      };
+      const auto key = ij2key(i,j);
+      const auto it = E2V.find(key);
+      if(it == E2V.end())
+      {
+        // new vertex
+        if(n==V.rows()){ V.conservativeResize(V.rows()*2+1,V.cols()); }
+        v = n;
+        E2V[key] = v;
+        n++;
+        return true;
+      }else
+      {
+        v = it->second;
+        return false;
+      }
+
     };
-    const auto key = ij2key(i,j);
-    const auto it = E2V.find(key);
-    int v = -1;
-    if(it == E2V.end())
-    {
-      // new vertex
-      if(n==V.rows()){ V.conservativeResize(V.rows()*2+1,V.cols()); }
-      V.row(n) = GV.row(i) + t*(GV.row(j) - GV.row(i));
-      v = n;
-      E2V[key] = v;
-      n++;
-    }else
-    {
-      v = it->second;
-    }
-    return v;
-  };
 
-    int c_flags = 0;
-    for(int c = 0; c < 8; c++)
-    {
-      if(cS(c) > isovalue){ c_flags |= 1<<c; }
-    }
-    //Find which edges are intersected by the surface
-    int e_flags = aiCubeEdgeFlags[c_flags];
-    //If the cube is entirely inside or outside of the surface, then there will be no intersections
-    if(e_flags == 0) { return; }
-    //Find the point of intersection of the surface with each edge
-    //Then find the normal to the surface at those points
-    Eigen::Matrix<Index,12,1> edge_vertices;
-    for(int e = 0; e < 12; e++)
-    {
+  int c_flags = 0;
+  for(int c = 0; c < 8; c++)
+  {
+    if(cS(c) > isovalue){ c_flags |= 1<<c; }
+  }
+  //Find which edges are intersected by the surface
+  int e_flags = aiCubeEdgeFlags[c_flags];
+  //If the cube is entirely inside or outside of the surface, then there will be no intersections
+  if(e_flags == 0) { return; }
+  //Find the point of intersection of the surface with each edge
+  //Then find the normal to the surface at those points
+  Eigen::Matrix<Index,12,1> edge_vertices;
+  for(int e = 0; e < 12; e++)
+  {
 #ifndef NDEBUG
-      edge_vertices[e] = -1;
+    edge_vertices[e] = -1;
 #endif
-      //if there is an intersection on this edge
-      if(e_flags & (1<<e))
+    //if there is an intersection on this edge
+    if(e_flags & (1<<e))
+    {
+      // record global index into local table
+      const Index & i = cI(a2eConnection[e][0]);
+      const Index & j = cI(a2eConnection[e][1]);
+      Index & v = edge_vertices[e];
+      if(ij2vertex(i,j,v))
       {
         // find crossing point assuming linear interpolation along edges
         const Scalar & a = cS(a2eConnection[e][0]);
         const Scalar & b = cS(a2eConnection[e][1]);
+        // This t is being recomputed each time an edge is seen.
         Scalar t;
-        {
-          const Scalar delta = b-a;
-          if(delta == 0) { t = 0.5; }
-          t = (isovalue - a)/delta;
-        };
-        // record global index into local table
-        edge_vertices[e] = 
-          ij2vertex(cI(a2eConnection[e][0]),cI(a2eConnection[e][1]),t);
-        assert(edge_vertices[e] >= 0);
-        assert(edge_vertices[e] < n);
+        const Scalar delta = b-a;
+        if(delta == 0) { t = 0.5; }
+        t = (isovalue - a)/delta;
+        V.row(v) = GV.row(i) + t*(GV.row(j) - GV.row(i));
       }
+      assert(v >= 0);
+      assert(v < n);
     }
-    // Insert the triangles that were found.  There can be up to five per cube
-    for(int f = 0; f < 5; f++)
-    {
-      if(a2fConnectionTable[c_flags][3*f] < 0) break;
-      if(m==F.rows()){ F.conservativeResize(F.rows()*2+1,F.cols()); }
-      assert(edge_vertices[a2fConnectionTable[c_flags][3*f+0]]>=0);
-      assert(edge_vertices[a2fConnectionTable[c_flags][3*f+1]]>=0);
-      assert(edge_vertices[a2fConnectionTable[c_flags][3*f+2]]>=0);
-      F.row(m) <<
-        edge_vertices[a2fConnectionTable[c_flags][3*f+0]],
-        edge_vertices[a2fConnectionTable[c_flags][3*f+1]],
-        edge_vertices[a2fConnectionTable[c_flags][3*f+2]];
+  }
+  // Insert the triangles that were found.  There can be up to five per cube
+  for(int f = 0; f < 5; f++)
+  {
+    if(a2fConnectionTable[c_flags][3*f] < 0) break;
+    if(m==F.rows()){ F.conservativeResize(F.rows()*2+1,F.cols()); }
+    assert(edge_vertices[a2fConnectionTable[c_flags][3*f+0]]>=0);
+    assert(edge_vertices[a2fConnectionTable[c_flags][3*f+1]]>=0);
+    assert(edge_vertices[a2fConnectionTable[c_flags][3*f+2]]>=0);
+    F.row(m) <<
+      edge_vertices[a2fConnectionTable[c_flags][3*f+0]],
+      edge_vertices[a2fConnectionTable[c_flags][3*f+1]],
+      edge_vertices[a2fConnectionTable[c_flags][3*f+2]];
       m++;
-    }
+  }
 }
 
 

include/igl/marching_cubes.cpp

@@ -25,6 +25,26 @@ IGL_INLINE void igl::marching_cubes(
     const typename DerivedS::Scalar isovalue,
     Eigen::PlainObjectBase<DerivedV> &V,
     Eigen::PlainObjectBase<DerivedF> &F)
+{
+  std::unordered_map<std::int64_t,int> E2V;
+  return marching_cubes(S,GV,nx,ny,nz,isovalue,V,F,E2V);
+}
+
+template <
+  typename DerivedS, 
+  typename DerivedGV, 
+  typename DerivedV, 
+  typename DerivedF>
+IGL_INLINE void igl::marching_cubes(
+  const Eigen::MatrixBase<DerivedS> & S,
+  const Eigen::MatrixBase<DerivedGV> & GV,
+  const unsigned nx,
+  const unsigned ny,
+  const unsigned nz,
+  const typename DerivedS::Scalar isovalue,
+  Eigen::PlainObjectBase<DerivedV> &V,
+  Eigen::PlainObjectBase<DerivedF> &F,
+  std::unordered_map<std::int64_t,int> &E2V)
 {
   typedef typename DerivedS::Scalar Scalar;
   typedef unsigned Index;
@@ -32,7 +52,6 @@ IGL_INLINE void igl::marching_cubes(
   const unsigned ioffset[8] = {0,1,1+nx,nx,nx*ny,1+nx*ny,1+nx+nx*ny,nx+nx*ny};
 
 
-  std::unordered_map<std::int64_t,int> E2V;
   V.resize(std::pow(nx*ny*nz,2./3.),3);
   F.resize(std::pow(nx*ny*nz,2./3.),3);
   Index n = 0;

include/igl/marching_cubes.h

@@ -8,6 +8,7 @@
 #ifndef IGL_MARCHING_CUBES_H
 #define IGL_MARCHING_CUBES_H
 #include "igl_inline.h"
+#include <unordered_map>
 
 #include <Eigen/Core>
 namespace igl
@@ -39,6 +40,27 @@ namespace igl
     const typename DerivedS::Scalar isovalue,
     Eigen::PlainObjectBase<DerivedV> &V,
     Eigen::PlainObjectBase<DerivedF> &F);
+  /// \overload 
+  /// 
+  /// \brief Return edge-to-vertex map which can be used to implement
+  /// batched root finding by caller (see 909_BatchMarchingCubes)
+  ///
+  /// @param[out] E2V  map from edge key to index into rows of V
+  template <
+    typename DerivedS, 
+    typename DerivedGV, 
+    typename DerivedV, 
+    typename DerivedF>
+  IGL_INLINE void marching_cubes(
+    const Eigen::MatrixBase<DerivedS> & S,
+    const Eigen::MatrixBase<DerivedGV> & GV,
+    const unsigned nx,
+    const unsigned ny,
+    const unsigned nz,
+    const typename DerivedS::Scalar isovalue,
+    Eigen::PlainObjectBase<DerivedV> &V,
+    Eigen::PlainObjectBase<DerivedF> &F,
+    std::unordered_map<std::int64_t,int> &E2V);
   /// \overload
   ///
   /// \brief Sparse voxel version

tutorial/909_BatchMarchingCubes/main.cpp

@@ -0,0 +1,126 @@
+#include <igl/marching_cubes.h>
+#include <igl/signed_distance.h>
+#include <igl/read_triangle_mesh.h>
+#include <igl/voxel_grid.h>
+#include <igl/opengl/glfw/Viewer.h>
+#include <Eigen/Core>
+#include <iostream>
+
+
+int main(int argc, char * argv[])
+{
+  using namespace Eigen;
+  using namespace std;
+  using namespace igl;
+  MatrixXi F;
+  MatrixXd V;
+  // Read in inputs as double precision floating point meshes
+  read_triangle_mesh(
+      TUTORIAL_SHARED_PATH "/armadillo.obj",V,F);
+  cout<<"Creating grid..."<<endl;
+  // number of vertices on the largest side
+  const int s = 100;
+  // create grid
+  MatrixXd GV;
+  Eigen::RowVector3i res;
+  igl::voxel_grid(V,0,s,1,GV,res);
+ 
+  // compute values
+  cout<<"Computing distances..."<<endl;
+
+  // Batch based function for evaluating implicit
+  const auto batch_implicit = [&V,&F](const MatrixXd & Q)->VectorXd
+  {
+    VectorXd S;
+    {
+      VectorXi I;
+      MatrixXd C,N;
+      signed_distance(Q,V,F,SIGNED_DISTANCE_TYPE_FAST_WINDING_NUMBER,S,I,C,N);
+      // Extremely flatten out near zero
+      S = S.array().sign() * S.array().abs().exp();
+    }
+    return S;
+  };
+
+  VectorXd S = batch_implicit(GV);
+  cout<<"Marching cubes..."<<endl;
+  MatrixXd SV;
+  MatrixXi SF;
+  igl::marching_cubes(S,GV,res(0),res(1),res(2),0,SV,SF);
+
+  std::unordered_map<std::int64_t,int> E2V;
+  igl::marching_cubes(S,GV,res(0),res(1),res(2),0,SV,SF,E2V);
+
+  // Initialize min and max for root finding bisection
+  assert(E2V.size() == SV.rows());
+  Eigen::MatrixXd T(SV.rows(),2);
+  Eigen::VectorXi I(SV.rows());
+
+  // Precompute slices for end points
+  Eigen::MatrixXd GVi(SV.rows(),3);
+  Eigen::MatrixXd GVj(SV.rows(),3);
+
+  // This is only used for the assertion below
+  const auto ij2key = [](std::int32_t i,std::int32_t j)
+  {
+    if(i>j){ std::swap(i,j); }
+    std::int64_t ret = 0;
+    ret |= i;
+    ret |= static_cast<std::int64_t>(j) << 32;
+    return ret;
+  };
+  const auto key2ij = [](const std::int64_t & key, std::int32_t & i, std::int32_t & j)
+  {
+    i = key & 0xFFFFFFFF;
+    j = key >> 32;
+  };
+  for (const auto& e2v: E2V)
+  {
+    const std::int64_t key = e2v.first;
+    const std::int32_t v = e2v.second;
+    std::int32_t i,j;
+    key2ij(key,i,j);
+    const std::int64_t key0 = ij2key(i,j);
+    assert(key0 != key);
+    T.row(v) << 0,1;
+    // (i,j) is ordered so that i<j, but let's order so that S(i)<S(j)
+    if(S(i)>S(j)) { std::swap(i,j); }
+    GVi.row(v) = GV.row(i);
+    GVj.row(v) = GV.row(j);
+  }
+
+  const auto root_find_iteration = [&SV,&I,&T,&GVi,&GVj,&batch_implicit]()
+  {
+    Eigen::VectorXd T_mid = (T.col(0)+T.col(1))/2;
+    // Use this midpoint guess for the current visualization
+    SV = GVi.array().colwise() * (1.0 - T_mid.array()) + GVj.array().colwise() * T_mid.array();
+    // Compute values at midpoints
+    VectorXd S_mid = batch_implicit(SV);
+    // Update bounds
+    T.col(1) = (S_mid.array() >  0).select(T_mid, T.col(1));
+    T.col(0) = (S_mid.array() <= 0).select(T_mid, T.col(0));
+  };
+
+  cout<<R"(Usage:
+' '  Conduct a bisection iteration
+)";
+  igl::opengl::glfw::Viewer viewer;
+  viewer.data().set_mesh(SV,SF);
+  viewer.data().show_lines = false;
+  viewer.callback_key_down =
+    [&](igl::opengl::glfw::Viewer & viewer, unsigned char key, int mod)->bool
+    {
+      switch(key)
+      {
+        default:
+          return false;
+        case ' ':
+          root_find_iteration();
+          viewer.data().set_vertices(SV);
+          viewer.data().compute_normals();
+          break;
+      }
+      return true;
+    };
+  viewer.launch();
+}

tutorial/CMakeLists.txt

@@ -128,4 +128,5 @@ if(LIBIGL_TUTORIALS_CHAPTER9)
     igl_add_tutorial(906_TrimWithSolid igl::glfw igl_copyleft::cgal)
     igl_add_tutorial(907_DynamicAABB igl::glfw)
     igl_add_tutorial(908_IntersectionBlockingDecimation igl::glfw igl::predicates)
+    igl_add_tutorial(909_BatchMarchingCubes igl::glfw)
 endif()