version with precomputed grid values

include/igl/dual_contouring.cpp

@@ -325,18 +325,19 @@ namespace igl
         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)
+        igl::parallel_for(nx,[&](const int x)
         {
-          const int x = xm1+1;
-          for(int y = 1;y<ny;y++)
+          for(int y = 0;y<ny;y++)
           {
-            for(int z = 1;z<nz;z++)
+            for(int z = 0;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++)
               {
+                // off-by-one boundary cases
+                if(((o==0)&&x==0)||((o==1)&&y==0)||((o==2)&&z==0)){ continue;}
                 single_edge(x,y,z,o,e0,f0);
               }
             }
@@ -344,6 +345,56 @@ namespace igl
         },10ul);
         dual_vertex_positions();
       }
+      template <typename DerivedGf, typename DerivedGV>
+      void dense(
+        const Eigen::MatrixBase<DerivedGf> & Gf,
+        const Eigen::MatrixBase<DerivedGV> & GV,
+        const int nx,
+        const int ny,
+        const int nz)
+      {
+        min_corner = GV.colwise().minCoeff();
+        const RowVector3S max_corner = GV.colwise().maxCoeff();
+        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);
+
+        const auto xyz2i = [&nx,&ny,&nz]
+          (const int & x, const int & y, const int & z)->Eigen::Index
+        {
+          return x+nx*(y+ny*(z));
+        };
+
+        // loop over grid
+        //for(int z = 0;z<nz;z++)
+        igl::parallel_for(nz,[&](const int z)
+        {
+          for(int y = 0;y<ny;y++)
+          {
+            for(int x = 0;x<nx;x++)
+            {
+              //const Scalar f0 = f(e0);
+              const Eigen::Index k0 = xyz2i(x,y,z);
+              const RowVector3S e0 = GV.row(k0);
+              const Scalar f0 = Gf(k0);
+              // we'll consider the edges going "back" from this vertex
+              for(int o = 0;o<3;o++)
+              {
+                Eigen::RowVector3i jc(x,y,z);
+                jc(o) -= 1;
+                if(jc(o)<0) { continue;} // off-by-one boundary cases
+                const int k1 = xyz2i(jc(0),jc(1),jc(2));
+                const RowVector3S e1 = GV.row(k1);
+                const Scalar f1 = Gf(k1);
+                single_edge(x,y,z,o,e0,f0,e1,f1);
+              }
+            }
+          }
+        },10ul);
+        dual_vertex_positions();
+      }
       void sparse(
         const RowVector3S & _step,
         const Eigen::Matrix<Scalar,Eigen::Dynamic,1> & Gf,
@@ -438,6 +489,35 @@ IGL_INLINE void igl::dual_contouring(
   Q = DC.Q.template cast<typename DerivedQ::Scalar>();
 }
 
+template <
+  typename DerivedGf,
+  typename DerivedGV,
+  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::MatrixBase<DerivedGf> & Gf,
+  const Eigen::MatrixBase<DerivedGV> & GV,
+  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(Gf,GV,nx,ny,nz);
+  V = DC.V;
+  Q = DC.Q.template cast<typename DerivedQ::Scalar>();
+}
+
 template <
   typename DerivedGf,
   typename DerivedGV,
@@ -469,7 +549,10 @@ IGL_INLINE void igl::dual_contouring(
 #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> >&);
+template void igl::dual_contouring<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 3, 1, -1, 3>, 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::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -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<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::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > 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> >&);
+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> >&);
+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, -1, 0, -1, -1>, 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, -1, 0, -1, -1> > const&, bool, bool, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 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

@@ -47,6 +47,32 @@ namespace igl
     const bool root_finding,
     Eigen::PlainObjectBase<DerivedV> & V,
     Eigen::PlainObjectBase<DerivedQ> & Q);
+  // Inputs:
+  //   Gf  nx*ny*nz list of function values so that Gf(k) = f(GV.row(k)) (only
+  //   needs to be accurate near f=0 and correct sign elsewhere)
+  //   GV  nx*ny*nz list of grid positions so that the x,y,z grid position is at
+  //     GV.row(x+nx*(y+z*ny))
+  template <
+    typename DerivedGf,
+    typename DerivedGV,
+    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::MatrixBase<DerivedGf> & Gf,
+    const Eigen::MatrixBase<DerivedGV> & GV,
+    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