dual contouring (dense/sparse); tutorial

include/igl/dual_contouring.cpp

@@ -0,0 +1,475 @@
+#include "dual_contouring.h"
+#include "quadprog.h"
+#include "parallel_for.h"
+#include <thread>
+#include <vector>
+#include <unordered_map>
+#include <iostream>
+
+namespace igl
+{
+  // These classes not intended to be used directly
+  class Hash
+  {
+    public:
+      // https://stackoverflow.com/a/26348708/148668
+      uint64_t operator()(const std::tuple<int, int, int> & key) const
+      {
+        // Check that conversion is safe. Could use int16_t directly everywhere
+        // below but it's an uncommon type to expose and grid indices should
+        // never be more than 2¹⁵-1 in the first place.
+        assert( std::get<0>(key) == (int)(int16_t)std::get<0>(key));
+        assert( std::get<1>(key) == (int)(int16_t)std::get<1>(key));
+        assert( std::get<2>(key) == (int)(int16_t)std::get<2>(key));
+        uint64_t result = uint16_t(std::get<0>(key));
+        result = (result << 16) + uint16_t(std::get<1>(key));
+        result = (result << 16) + uint16_t(std::get<2>(key));
+        return result;
+      };
+  };
+  template <typename Scalar>
+  class DualContouring
+  {
+    // Types
+    public:
+      typedef Eigen::Matrix<Scalar,1,3> RowVector3S;
+      typedef Eigen::Matrix<Scalar,1,4> RowVector4S;
+      typedef Eigen::Matrix<Scalar,4,4> Matrix4S;
+      typedef Eigen::Matrix<Scalar,3,3> Matrix3S;
+      typedef Eigen::Matrix<Scalar,3,1> Vector3S;
+      typedef std::tuple<int, int, int> KeyTriplet;
+    // helper function
+    public:
+    // Working variables
+    public: 
+      // f(x) returns >0 outside, <0 inside, and =0 on the surface
+      std::function<Scalar(const RowVector3S &)> f;
+      // f_grad(x) returns (df/dx)/‖df/dx‖ (normalization only important when
+      // f(x) = 0).
+      std::function<RowVector3S(const RowVector3S &)> f_grad;
+      bool constrained;
+      bool triangles;
+      bool root_finding;
+      RowVector3S min_corner;
+      RowVector3S step;
+      Eigen::Matrix<Scalar,Eigen::Dynamic,3> V;
+      // Running number of vertices added during contouring
+      typename decltype(V)::Index n;
+      std::unordered_map< KeyTriplet, typename decltype(V)::Index, Hash > C2V;
+      std::vector<RowVector3S,Eigen::aligned_allocator<RowVector3S>> vV;
+      std::vector<Eigen::RowVector3i,Eigen::aligned_allocator<Eigen::RowVector3i>> vI;
+      std::vector<Matrix4S,Eigen::aligned_allocator<Matrix4S>> vH;
+      std::vector<int> vcount;
+      Eigen::Matrix<Eigen::Index,Eigen::Dynamic,Eigen::Dynamic> Q;
+      typename decltype(Q)::Index m;
+      std::mutex Qmut;
+      std::mutex Vmut;
+    public:
+      DualContouring(
+        const std::function<Scalar(const RowVector3S &)> & _f,
+        const std::function<RowVector3S(const RowVector3S &)> & _f_grad,
+        const bool _constrained = false,
+        const bool _triangles = false,
+        const bool _root_finding = true):
+        f(_f),
+        f_grad(_f_grad),
+        constrained(_constrained),
+        triangles(_triangles),
+        root_finding(_root_finding),
+        n(0),
+        C2V(0),
+        vV(),vI(),vH(),vcount(),
+        m(0)
+      { }
+      // Side effects: new entry in vV,vI,vH,vcount, increment n
+      // Returns index of new vertex
+      typename decltype(V)::Index new_vertex()
+      {
+        const auto v = n;
+        n++;
+        vcount.resize(n);
+        vV.resize(n);
+        vI.resize(n);
+        vH.resize(n);
+        vcount[v] = 0;
+        vV[v].setZero();
+        vH[v].setZero();
+        return v;
+      };
+      // Inputs:
+      //   kc  3-long vector of {x,y,z} index of primal grid **cell**
+      // Returns index to corresponding dual vertex
+      // Side effects: if vertex for this cell does not yet exist, creates it
+      typename decltype(V)::Index sub2dual(const Eigen::RowVector3i & kc)
+      {
+        const KeyTriplet key = {kc(0),kc(1),kc(2)};
+        const auto it = C2V.find(key);
+        typename decltype(V)::Index v = -1;
+        if(it == C2V.end())
+        {
+          v = new_vertex();
+          C2V[key] = v;
+          vI[v] = kc;
+        }else
+        {
+          v = it->second;
+        }
+        return v;
+      };
+      RowVector3S primal(const Eigen::RowVector3i & ic) const
+      {
+        return min_corner + (ic.cast<Scalar>().array() * step.array()).matrix();
+      }
+      Eigen::RowVector3i inverse_primal(const RowVector3S & x) const
+      {
+        // x = min_corner + (ic.cast<Scalar>().array() * step.array()).matrix();
+        // x-min_corner = (ic.cast<Scalar>().array() * step.array())
+        // (x-min_corner).array() / step.array()  = ic.cast<Scalar>().array()
+        // ((x-min_corner).array() / step.array()).round()  = ic
+        return 
+          ((x-min_corner).array()/step.array()).round().template cast<int>();
+      }
+      // Inputs:
+      //   x  x-index of vertex on primal grid
+      //   y  y-index of vertex on primal grid
+      //   z  z-index of vertex on primal grid
+      //   o  which edge are we looking back on? o=0->x,o=1->y,o=2->z
+      // Side effects: may insert new vertices into vV,vI,vH,vcount, new faces
+      //   into Q
+      bool single_edge(const int & x, const int & y, const int & z, const int & o)
+      {
+        const RowVector3S e0 = primal(Eigen::RowVector3i(x,y,z));
+        const Scalar f0 = f(e0);
+        return single_edge(x,y,z,o,e0,f0);
+      }
+      bool single_edge(
+        const int & x, 
+        const int & y, 
+        const int & z, 
+        const int & o,
+        const RowVector3S & e0,
+        const Scalar & f0)
+      {
+        //e1 computed here needs to precisely agree with e0 when called with
+        //correspond x,y,z. So, don't do this:
+        //Eigen::RowVector3d e1 = e0;
+        //e1(o) -= step(o);
+        Eigen::RowVector3i jc(x,y,z);
+        jc(o) -= 1;
+        const RowVector3S e1 = primal(jc);
+        const Scalar f1 = f(e1);
+        return single_edge(x,y,z,o,e0,f0,e1,f1);
+      }
+      bool single_edge(
+        const int & x, 
+        const int & y, 
+        const int & z, 
+        const int & o,
+        const RowVector3S & e0,
+        const Scalar & f0,
+        const RowVector3S & e1,
+        const Scalar & f1)
+      {
+        const Scalar isovalue = 0;
+        if((f0>isovalue) == (f1>isovalue)) { return false; }
+        // Position of crossing point along edge
+        RowVector3S p;
+        Scalar t = -1;
+        if(root_finding)
+        {
+          Scalar tl = 0;
+          bool gl = f0>0;
+          Scalar tu = 1;
+          bool gu = f1>0;
+          assert(gu ^ gl);
+          int riter = 0;
+          const int max_riter = 7;
+          while(true)
+          {
+            t = 0.5*(tu + tl);
+            p = e0+t*(e1-e0);
+            riter++;
+            if(riter > max_riter) { break;}
+            const Scalar ft = f(p);
+            if( (ft>0) == gu) { tu = t; }
+            else if( (ft>0) == gl){ tl = t; }
+            else { break; }
+          }
+        }else
+        {
+          // inverse lerp
+          const Scalar delta = f1-f0;
+          if(delta == 0) { t = 0.5; }
+          t = (isovalue - f0)/delta;
+          p = e0+t*(e1-e0);
+        }
+        // insert vertex at this point to triangulate quad face
+        const typename decltype(V)::Index ev = triangles ? new_vertex() : -1;
+        if(triangles)
+        {
+          const std::lock_guard<std::mutex> lock(Vmut);
+          vV[ev] = p;
+          vcount[ev] = 1;
+          vI[ev] = Eigen::RowVector3i(-1,-1,-1);
+        }
+        // edge normal from function handle (could use grid finite
+        // differences/interpolation gradients)
+        const RowVector3S dfdx = f_grad(p);
+        // homogenous plane equation
+        const RowVector4S P = (RowVector4S()<<dfdx,-dfdx.dot(p)).finished();
+        // quadric contribution
+        const Matrix4S H = P.transpose() * P;
+        // Build quad face from dual vertices of 4 cells around this edge
+        Eigen::RowVector4i face;
+        int k = 0;
+        for(int i = -1;i<=0;i++)
+        {
+          for(int j = -1;j<=0;j++)
+          {
+            Eigen::RowVector3i kc(x,y,z);
+            kc(o)--;
+            kc((o+1)%3)+=i;
+            kc((o+2)%3)+=j;
+            const std::lock_guard<std::mutex> lock(Vmut);
+            const typename decltype(V)::Index v = sub2dual(kc);
+            vV[v] += p;
+            vcount[v]++;
+            vH[v] += H;
+            face(k++) = v;
+          }
+        }
+        {
+          const std::lock_guard<std::mutex> lock(Qmut);
+          if(triangles)
+          {
+            if(m+4 >= Q.rows()){ Q.conservativeResize(2*m+4,Q.cols()); }
+            if(f0>f1)
+            {
+              Q.row(m+0)<<      ev,face(3),face(1)        ;
+              Q.row(m+1)<<              ev,face(1),face(0);
+              Q.row(m+2)<< face(2),             ev,face(0);
+              Q.row(m+3)<< face(2),face(3),             ev;
+            }else
+            {
+              Q.row(m+0)<<      ev,face(1),face(3)        ;
+              Q.row(m+1)<<              ev,face(3),face(2);
+              Q.row(m+2)<< face(0),             ev,face(2);
+              Q.row(m+3)<< face(0),face(1),             ev;
+            }
+            m+=4;
+          }else
+          {
+            if(m+1 >= Q.rows()){ Q.conservativeResize(2*m+1,Q.cols()); }
+            if(f0>f1)
+            {
+              Q.row(m)<< face(2),face(3),face(1),face(0);
+            }else
+            {
+              Q.row(m)<< face(0),face(1),face(3),face(2);
+            }
+            m++;
+          }
+        }
+        return true;
+      }
+      // Side effects: Q resized to fit m, V constructed to fit n and
+      // reconstruct data in vH,vI,vV,vcount
+      void dual_vertex_positions()
+      {
+        Q.conservativeResize(m,Q.cols());
+        V.resize(n,3);
+        //for(int v = 0;v<n;v++)
+        igl::parallel_for(n,[&](const Eigen::Index v)
+        {
+          RowVector3S mid = vV[v] / Scalar(vcount[v]);
+          if(triangles && vI[v](0)<0 ){ V.row(v) = mid; return; }
+          const Scalar w = 1e-2*(0.01+vcount[v]);
+          Matrix3S A = vH[v].block(0,0,3,3) + w*Matrix3S::Identity();
+          RowVector3S b = -vH[v].block(3,0,1,3) + w*mid;
+          // Replace with solver
+          //RowVector3S p = b * A.inverse();
+          //
+          // min_p  ½ pᵀ A p - pᵀb
+          //  
+          // let p = p₀ + x
+          //
+          //     min    ½ (p₀ + x )ᵀ A (p₀ + x ) - (p₀ + x )ᵀb
+          // step≥x≥0
+          const RowVector3S p0 = 
+            min_corner + ((vI[v].template cast<Scalar>().array()) * step.array()).matrix();
+          const RowVector3S x = 
+            constrained ?
+            igl::quadprog<Scalar,3>(A,(p0*A-b).transpose(),Vector3S(0,0,0),step.transpose()) :
+            Eigen::LLT<Matrix3S>(A).solve(-(p0*A-b).transpose());
+          V.row(v) = p0+x;
+        },1000ul);
+      }
+      // Inputs:
+      //   _min_corner  minimum (bottomLeftBack) corner of primal grid
+      //   max_corner  maximum (topRightFront) corner of primal grid
+      //   nx  number of primal grid vertices along x-axis
+      //   ny  number of primal grid vertices along y-ayis
+      //   nz  number of primal grid vertices along z-azis
+      // Side effects: prepares vV,vI,vH,vcount, Q for vertex_positions()
+      void dense(
+        const RowVector3S & _min_corner,
+        const RowVector3S & max_corner,
+        const int nx,
+        const int ny,
+        const int nz)
+      {
+        min_corner = _min_corner;
+        step =
+          (max_corner-min_corner).array()/(RowVector3S(nx,ny,nz).array()-1);
+        // Should do some reasonable reserves for C2V,vV,vI,vH,vcount
+        Q.resize(std::pow(nx*ny*nz,2./3.),triangles?3:4);
+        // loop over grid
+        //for(int x = 1;x<nx;x++)
+        igl::parallel_for(nx-1,[&](const int xm1)
+        {
+          const int x = xm1+1;
+          for(int y = 1;y<ny;y++)
+          {
+            for(int z = 1;z<nz;z++)
+            {
+              const RowVector3S e0 = primal(Eigen::RowVector3i(x,y,z));
+              const Scalar f0 = f(e0);
+              // we'll consider the edges going "back" from this vertex
+              for(int o = 0;o<3;o++)
+              {
+                single_edge(x,y,z,o,e0,f0);
+              }
+            }
+          }
+        },10ul);
+        dual_vertex_positions();
+      }
+      void sparse(
+        const RowVector3S & _step,
+        const Eigen::Matrix<Scalar,Eigen::Dynamic,1> & Gf,
+        const Eigen::Matrix<Scalar,Eigen::Dynamic,3> & GV,
+        const Eigen::Matrix<int,Eigen::Dynamic,2> & GI)
+      {
+        step = _step;
+        Q.resize((triangles?4:1)*GI.rows(),triangles?3:4);
+        // in perfect world doesn't matter where min_corner is so long as it is
+        // _on_ the grid. For models very far from origin, centering grid near
+        // model avoids possible rounding error in hash()/inverse_primal()
+        // [still very unlikely, but let's be safe]
+        min_corner = GV.colwise().minCoeff();
+        // igl::parallel_for here made things worse. Probably need to do proper
+        // map-reduce rather than locks on mutexes.
+        for(Eigen::Index i = 0;i<GI.rows();i++)
+        {
+          RowVector3S e0 = GV.row(GI(i,0));
+          RowVector3S e1 = GV.row(GI(i,1));
+          Eigen::RowVector3i ic0 = inverse_primal(e0);
+          Eigen::RowVector3i ic1 = inverse_primal(e1);
+#ifndef NDEBUG
+          RowVector3S p0 = primal(ic0);
+          RowVector3S p1 = primal(ic1);
+          assert( (p0-e0).norm() < 1e-10);
+          assert( (p1-e1).norm() < 1e-10);
+#endif
+          Scalar f0 = Gf(GI(i,0)); //f(e0);
+          Scalar f1 = Gf(GI(i,1)); //f(e1);
+          // should differ in just one coordinate. Find that coordinate.
+          int o = -1;
+          for(int j = 0;j<3;j++)
+          {
+            if(ic0(j) == ic1(j)){ continue;}
+            if(ic0(j) - ic1(j) == 1)
+            { 
+              assert(o == -1 && "Edges should differ in just one coordinate"); 
+              o = j; 
+              continue; // rather than break so assertions fire
+            }
+            if(ic1(j) - ic0(j) == 1)
+            { 
+              assert(o == -1 && "Edges should differ in just one coordinate"); 
+              std::swap(e0,e1);
+              std::swap(f0,f1);
+              std::swap(ic0,ic1);
+              o = j; 
+              continue; // rather than break so assertions fire
+            } else
+            {
+              assert(false && "Edges should differ in just one coordinate"); 
+            }
+          }
+          assert(o>=0 && "Edges should differ in just one coordinate");
+          // i0 is the larger subscript location and ic1 is backward in the o
+          // direction.
+          for(int j = 0;j<3;j++){ assert(ic0(j) == ic1(j)+(o==j)); } 
+          const int x = ic0(0);
+          const int y = ic0(1);
+          const int z = ic0(2);
+          single_edge(x,y,z,o,e0,f0,e1,f1);
+        }
+        dual_vertex_positions();
+      }
+  };
+}
+
+template <
+  typename DerivedV,
+  typename DerivedQ>
+IGL_INLINE void igl::dual_contouring(
+  const std::function<
+    typename DerivedV::Scalar(const Eigen::Matrix<typename DerivedV::Scalar,1,3> &)> & f,
+  const std::function<
+    Eigen::Matrix<typename DerivedV::Scalar,1,3>(
+      const Eigen::Matrix<typename DerivedV::Scalar,1,3> &)> & f_grad,
+  const Eigen::Matrix<typename DerivedV::Scalar,1,3> & min_corner,
+  const Eigen::Matrix<typename DerivedV::Scalar,1,3> & max_corner,
+  const int nx,
+  const int ny,
+  const int nz,
+  const bool constrained,
+  const bool triangles,
+  const bool root_finding,
+  Eigen::PlainObjectBase<DerivedV> & V,
+  Eigen::PlainObjectBase<DerivedQ> & Q)
+{
+  typedef typename DerivedV::Scalar Scalar;
+  DualContouring<Scalar> DC(f,f_grad,constrained,triangles,root_finding);
+  DC.dense(min_corner,max_corner,nx,ny,nz);
+  V = DC.V;
+  Q = DC.Q.template cast<typename DerivedQ::Scalar>();
+}
+
+template <
+  typename DerivedGf,
+  typename DerivedGV,
+  typename DerivedGI,
+  typename DerivedV,
+  typename DerivedQ>
+IGL_INLINE void igl::dual_contouring(
+  const std::function<typename DerivedV::Scalar(const Eigen::Matrix<typename DerivedV::Scalar,1,3> &)> & f,
+  const std::function<Eigen::Matrix<typename DerivedV::Scalar,1,3>(const Eigen::Matrix<typename DerivedV::Scalar,1,3> &)> & f_grad,
+  const Eigen::Matrix<typename DerivedV::Scalar,1,3> & step,
+  const Eigen::MatrixBase<DerivedGf> & Gf,
+  const Eigen::MatrixBase<DerivedGV> & GV,
+  const Eigen::MatrixBase<DerivedGI> & GI,
+  const bool constrained,
+  const bool triangles,
+  const bool root_finding,
+  Eigen::PlainObjectBase<DerivedV> & V,
+  Eigen::PlainObjectBase<DerivedQ> & Q)
+{
+  if(GI.rows() == 0){ return;}
+  assert(GI.cols() == 2);
+  typedef typename DerivedV::Scalar Scalar;
+  DualContouring<Scalar> DC(f,f_grad,constrained,triangles,root_finding);
+  DC.sparse(step,Gf,GV,GI);
+  V = DC.V;
+  Q = DC.Q.template cast<typename DerivedQ::Scalar>();
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::dual_contouring<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::function<Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar (Eigen::Matrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, 1, 3, 1, 1, 3> const&)> const&, std::function<Eigen::Matrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, 1, 3, 1, 1, 3> (Eigen::Matrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, 1, 3, 1, 1, 3> const&)> const&, Eigen::Matrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, 1, 3, 1, 1, 3> const&, Eigen::Matrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, 1, 3, 1, 1, 3> const&, int, int, int, bool, bool, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
+template void igl::dual_contouring<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::function<Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar (Eigen::Matrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, 1, 3, 1, 1, 3> const&)> const&, std::function<Eigen::Matrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, 1, 3, 1, 1, 3> (Eigen::Matrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, 1, 3, 1, 1, 3> const&)> const&, Eigen::Matrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, 1, 3, 1, 1, 3> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, bool, bool, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+#endif

include/igl/dual_contouring.h

@@ -0,0 +1,77 @@
+#ifndef IGL_DUAL_CONTOURING_H
+#define IGL_DUAL_CONTOURING_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <Eigen/Dense>
+#include <functional>
+namespace igl
+{
+  // Dual contouring to extract a pure quad mesh from differentiable implicit
+  // function using a dense grid.
+  //
+  // Inputs:
+  //   f  function returning >0 outside, <0 inside and =0 on the surface
+  //   f_grad  function returning ∇f/‖∇f‖
+  //   min_corner  position of primal grid vertex at minimum corner
+  //   max_corner  position of primal grid vertex at maximum corner
+  //   nx  number of vertices on x side of primal grid
+  //   ny  number of vertices on y side of primal grid
+  //   nz  number of vertices on z side of primal grid
+  //   constrained  whether to force dual vertices to lie strictly inside
+  //     corresponding primal cell (prevents self-intersections at cost of
+  //     surface quality; marginally slower)
+  //   triangles  whether to output four triangles instead of one quad per
+  //     crossing edge (quad mesh usually looks fine)
+  //   root_finding  whether to use root finding to identify crossing point on
+  //     each edge (improves quality a lot at cost of performance). If false,
+  //     use linear interpolation.
+  // Outputs:
+  //   V  #V by 3 list of outputs vertex positions
+  //   Q  #Q by 4 (or 3 if triangles=true) face indices into rows of V
+  template <
+    typename DerivedV,
+    typename DerivedQ>
+  IGL_INLINE void dual_contouring(
+    const std::function<
+      typename DerivedV::Scalar(const Eigen::Matrix<typename DerivedV::Scalar,1,3> &)> & f,
+    const std::function<
+      Eigen::Matrix<typename DerivedV::Scalar,1,3>(
+        const Eigen::Matrix<typename DerivedV::Scalar,1,3> &)> & f_grad,
+    const Eigen::Matrix<typename DerivedV::Scalar,1,3> & min_corner,
+    const Eigen::Matrix<typename DerivedV::Scalar,1,3> & max_corner,
+    const int nx,
+    const int ny,
+    const int nz,
+    const bool constrained,
+    const bool triangles,
+    const bool root_finding,
+    Eigen::PlainObjectBase<DerivedV> & V,
+    Eigen::PlainObjectBase<DerivedQ> & Q);
+  // Sparse voxel grid
+  //
+  // Gf  #GV list of corresponding f values. If using root finding then only the
+  //   sign needs to be correct.
+  template <
+    typename DerivedGf,
+    typename DerivedGV,
+    typename DerivedGI,
+    typename DerivedV,
+    typename DerivedQ>
+  IGL_INLINE void dual_contouring(
+    const std::function<typename DerivedV::Scalar(const Eigen::Matrix<typename DerivedV::Scalar,1,3> &)> & f,
+    const std::function<Eigen::Matrix<typename DerivedV::Scalar,1,3>(const Eigen::Matrix<typename DerivedV::Scalar,1,3> &)> & f_grad,
+    const Eigen::Matrix<typename DerivedV::Scalar,1,3> & step,
+    const Eigen::MatrixBase<DerivedGf> & Gf,
+    const Eigen::MatrixBase<DerivedGV> & GV,
+    const Eigen::MatrixBase<DerivedGI> & GI,
+    const bool constrained,
+    const bool triangles,
+    const bool root_finding,
+    Eigen::PlainObjectBase<DerivedV> & V,
+    Eigen::PlainObjectBase<DerivedQ> & Q);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "dual_contouring.cpp"
+#endif
+#endif

tutorial/715_MeshImplicitFunction/main.cpp

@@ -1,50 +1,242 @@
 #include <igl/copyleft/marching_cubes.h>
+#include <igl/unique_simplices.h>
+#include <igl/dual_contouring.h>
+#include <igl/get_seconds.h>
+#include <igl/grid.h>
+#include <igl/marching_cubes.h>
+#include <igl/per_corner_normals.h>
+#include <igl/parallel_for.h>
+#include <igl/per_corner_normals.h>
+#include <igl/per_face_normals.h>
+#include <igl/polygon_corners.h>
+#include <igl/slice.h>
 #include <igl/sparse_voxel_grid.h>
 #include <igl/opengl/glfw/Viewer.h>
 
-#include <Eigen/Core>
-#include <iostream>
-
-#include "tutorial_shared_path.h"
-
 int main(int argc, char * argv[])
 {
-  // An implicit function which is zero on the surface of a sphere centered at the origin with radius 1
-  // This function is negative inside the surface and positive outside the surface
-  std::function<double(const Eigen::RowVector3d&)> scalar_func = [](const Eigen::RowVector3d& pt) -> double {
-    return pt.norm() - 1.0;
+  const auto & tictoc = []()
+  {
+    static double t_start = igl::get_seconds();
+    double diff = igl::get_seconds()-t_start;
+    t_start += diff;
+    return diff;
+  };
+
+  // Create an interesting shape with sharp features using SDF CSG with spheres.
+  const auto & sphere = [](
+      const Eigen::RowVector3d & c,
+      const double r,
+      const Eigen::RowVector3d & x)->double
+  {
+    return (x-c).norm() - r;
+  };
+  const std::function<double(const Eigen::RowVector3d & x)> f = 
+    [&](const Eigen::RowVector3d & x)->double
+  {
+    return 
+      std::min(
+      std::min(std::max(
+      sphere(Eigen::RowVector3d(-0.2,0,-0.2),0.5,x),
+      -sphere(Eigen::RowVector3d(+0.2,0,0.2),0.5,x)),
+       sphere(Eigen::RowVector3d(-0.15,0,-0.15),0.3,x)
+      ),
+      std::max(
+      std::max(
+        sphere(Eigen::RowVector3d(-0.2,-0.5,-0.2),0.6,x),x(1)+0.45),-0.6-x(1))
+      );
   };
+  Eigen::RowVector3d p0(-0.2,0.5,-0.2);
+  assert(abs(f(p0)) < 1e-10 && "p0 should be on zero level-set");
 
-  // We know that the point (0, 0, 1) lies on the implicit surface
-  Eigen::RowVector3d p0(0., 0., 1.);
+  // Simple finite difference gradients
+  const auto & fd = [](
+    const std::function<double(const Eigen::RowVector3d&)> & f,
+    const Eigen::RowVector3d & x)
+  {
+    const double eps = 1e-10;
+    Eigen::RowVector3d g;
+    for(int c = 0;c<3;c++)
+    {
+      const Eigen::RowVector3d xp = x+eps*Eigen::RowVector3d(c==0,c==1,c==2);
+      const double fp = f(xp);
+      const Eigen::RowVector3d xn = x-eps*Eigen::RowVector3d(c==0,c==1,c==2);
+      const double fn = f(xn);
+      g(c) = (fp-fn)/(2*eps);
+    }
+    return g;
+  };
+  const auto & f_grad = [&fd,&f](const Eigen::RowVector3d & x)
+  {
+    return fd(f,x).normalized();
+  };
 
-  // Construct a sparse voxel grid whose cubes have edge length eps = 0.1.
-  // The cubes will form a thin shell around the implicit surface
-  const double eps = 0.1;
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi Q,F;
+  Eigen::MatrixXd mcV,mcN;
+  Eigen::MatrixXi mcF;
 
-  // CS will hold one scalar value at each cube vertex corresponding
-  // the value of the implicit at that vertex
-  Eigen::VectorXd CS;
+  // Grid parameters
+  const Eigen::RowVector3d min_corner(-2,-2,-2);
+  const Eigen::RowVector3d max_corner(+2,+2,+2);
+  const int s = 256;
+  int nx = s+1;
+  int ny = s+1;
+  int nz = s+1;
+  const Eigen::RowVector3d step =
+    (max_corner-min_corner).array()/(Eigen::RowVector3d(nx,ny,nz).array()-1);
+  // Sparse grid below assumes regular grid
+  assert((step(0) == step(1))&&(step(0) == step(2)));
 
-  // CV will hold the positions of the corners of the sparse voxel grid
-  Eigen::MatrixXd CV;
+  // Dual contouring parameters
+  bool constrained = false;
+  bool triangles = false;
+  bool root_finding = true;
+  for(int pass = 0;pass<2;pass++)
+  {
+    const bool sparse = pass == 1;
+    printf("Using %s grid..\n",sparse?"sparse":"dense");
+    if(sparse)
+    {
+      // igl::sparse_voxel_grid assumes (0,0,0) lies on the grid. But dense igl::grid
+      // below won't necessarily do that depending on nx,ny,nz.
+      tictoc();
+      Eigen::MatrixXd GV;
+      Eigen::VectorXd Gf;
+      Eigen::Matrix<int,Eigen::Dynamic,8> GI;
+      igl::sparse_voxel_grid(p0,f,step(0),16.*pow(step(0),-2.),Gf,GV,GI);
+      const auto t_Gf = tictoc();
+      printf("  %5f secs to populate sparse grid of %ld cells\n",t_Gf+tictoc(),GI.rows());
+      // Dual contouring requires list of sparse edges (not cells)
+      // extract _all_ edges from sparse_voxel_grid (conservative)
+      Eigen::Matrix<int,Eigen::Dynamic,2> GI2;
+      {
+        Eigen::Matrix<int,Eigen::Dynamic,2> all_GI2(GI.rows()*12,2);
+        all_GI2 << 
+          // front
+          GI.col(0),GI.col(1),
+          GI.col(1),GI.col(2),
+          GI.col(2),GI.col(3),
+          GI.col(3),GI.col(0),
+          // back
+          GI.col(4+0),GI.col(4+1),
+          GI.col(4+1),GI.col(4+2),
+          GI.col(4+2),GI.col(4+3),
+          GI.col(4+3),GI.col(4+0),
+          // sides
+          GI.col(0),GI.col(4+0), 
+          GI.col(1),GI.col(4+1), 
+          GI.col(2),GI.col(4+2),
+          GI.col(3),GI.col(4+3);
+        Eigen::VectorXi _1,_2;
+        igl::unique_simplices(all_GI2,GI2,_1,_2);
+      }
+      tictoc();
+      Eigen::RowVector3d step =
+        (max_corner-min_corner).array()/(Eigen::RowVector3d(nx,ny,nz).array()-1);
+      igl::dual_contouring(
+        f,f_grad,step,Gf,GV,GI2,constrained,triangles,root_finding,V,Q);
+      printf("  %5f secs dual contouring\n",t_Gf+tictoc());
+      // Could use igl::marching_cubes once
+      // https://github.com/libigl/libigl/pull/1687 is merged
+      tictoc();
+      igl::copyleft::marching_cubes(Gf,GV,GI,mcV, mcF);
+      printf("  %5f secs marching cubes\n",t_Gf+tictoc());
+    }else
+    {
 
-  // CI is a #cubes x 8 matrix of indices where each row contains the
-  // indices into CV of the 8 corners of a cube
-  Eigen::MatrixXi CI;
+      tictoc();
+      igl::dual_contouring(
+        f,f_grad,min_corner,max_corner,nx,ny,nz,constrained,triangles,root_finding,V,Q);
+      printf("  %5f secs dual contouring\n",tictoc());
+      // build and sample grid
+      tictoc();
+      Eigen::MatrixXd GV;
+      igl::grid(Eigen::RowVector3i(nx,ny,nz),GV);
+      Eigen::VectorXd Gf(GV.rows());
+      igl::parallel_for(GV.rows(),[&](const int i)
+      {
+        GV.row(i).array() *= (max_corner-min_corner).array();
+        GV.row(i) += min_corner;
+        Gf(i) = f(GV.row(i));
+      },1000ul);
+      const auto t_grid = tictoc();
+      igl::marching_cubes(Gf,GV,nx,ny,nz,0,mcV,mcF);
+      const auto t_mc = tictoc();
+      printf("  %5f secs (%5f + %5f) marching cubes\n",t_grid+t_mc,t_grid,t_mc);
+    }
+  }
 
-  // Construct the voxel grid, populating CS, CV, and CI
-  igl::sparse_voxel_grid(p0, scalar_func, eps, 1024 /*expected_number_of_cubes*/, CS, CV, CI);
+  // Crisp (as possible) rendering of resulting MC triangle mesh
+  igl::per_corner_normals(mcV,mcF,20,mcN);
+  // Crisp rendering of resulting DC quad mesh with edges
+  Eigen::MatrixXi E;
+  Eigen::MatrixXd VV,N,NN;
+  Eigen::VectorXi J;
+  Eigen::MatrixXi FF;
+  if(triangles)
+  {
+    VV = V;
+    FF = Q;
+    E.resize(Q.rows()*3,2);
+    E<<
+      Q.col(0), Q.col(1), 
+      Q.col(1), Q.col(2), 
+      Q.col(2), Q.col(0);
+  }else
+  {
+    Eigen::VectorXi I,C;
+    igl::polygon_corners(Q,I,C);
+    E.resize(Q.rows()*4,2);
+    E<<
+      Q.col(0), Q.col(1), 
+      Q.col(1), Q.col(2), 
+      Q.col(2), Q.col(3), 
+      Q.col(3), Q.col(0);
+    igl::per_face_normals(V,I,C,N,VV,FF,J);
+    igl::slice(N,J,1,NN);
+    igl::per_corner_normals(V,I,C,20,N,VV,FF,J,NN);
+  }
 
-  // Given the sparse voxel grid, use Marching Cubes to construct a triangle mesh of the surface
-  Eigen::MatrixXi F;
-  Eigen::MatrixXd V;
-  igl::copyleft::marching_cubes(CS, CV, CI, V, F);
-
-  // Draw the meshed implicit surface
-  igl::opengl::glfw::Viewer viewer;
-  viewer.data().clear();
-  viewer.data().set_mesh(V,F);
-  viewer.data().set_face_based(true);
-  viewer.launch();
+  igl::opengl::glfw::Viewer vr;
+  bool show_edges = true;
+  bool use_dc = true;
+  const auto update = [&]()
+  {
+    const bool was_face_based = vr.data().face_based ;
+    vr.data().clear();
+    if(use_dc)
+    {
+      vr.data().set_mesh(VV,FF);
+      vr.data().show_lines = false;
+      vr.data().set_normals(NN);
+      if(show_edges)
+      {
+        vr.data().clear_edges();
+        vr.data().set_edges(V,E,Eigen::RowVector3d(0,0,0));
+      }
+    }else
+    {
+      vr.data().set_mesh(mcV,mcF);
+      vr.data().set_normals(mcN);
+      vr.data().show_lines = show_edges;
+    }
+    vr.data().face_based = was_face_based;
+  };
+  update();
+  vr.data().face_based = true;
+  vr.callback_key_pressed = [&](decltype(vr) &,unsigned int key, int mod)
+  {
+    switch(key)
+    {
+      case ' ': use_dc=!use_dc; update();return true;
+      case 'L': case 'l': show_edges=!show_edges; update();return true;
+    }
+    return false;
+  };
+  std::cout<<R"(
+[space]  Toggle between dual contouring and marching cubes
+)";
+  vr.launch();
 }
+