Merge pull request #1716 from libigl/dual-contouring

Dual contouring

include/igl/copyleft/quadprog.cpp

@@ -1,5 +1,8 @@
 #include "quadprog.h"
+#include "../matlab_format.h"
 #include <vector>
+#include <iostream>
+#include <cstdio>
 /*
  FILE eiquadprog.hh
  
@@ -92,6 +95,22 @@ IGL_INLINE bool igl::copyleft::quadprog(
 {
   using namespace Eigen;
   typedef double Scalar;
+
+
+  const auto print_ivector= [](const char* name, const Eigen::MatrixXi & A, int n)
+  {
+    std::cout<<igl::matlab_format(A,name)<<std::endl;
+  };
+  const auto print_matrix = [](const char* name, const Eigen::MatrixXd & A, int n)
+  {
+    std::cout<<igl::matlab_format(A,name)<<std::endl;
+  };
+
+  const auto print_vector = [](const char* name, const Eigen::VectorXd & v, int n)
+  {
+    std::cout<<igl::matlab_format(v,name)<<std::endl;
+  };
+
   const auto distance = [](Scalar a, Scalar b)->Scalar
   {
   	Scalar a1, b1, t;

include/igl/cotmatrix.cpp

@@ -84,6 +84,7 @@ IGL_INLINE void igl::cotmatrix(
 #include "diag.h"
 #include "massmatrix.h"
 #include <Eigen/Geometry>
+#include <Eigen/QR>
 
 template <
   typename DerivedV, 

include/igl/decimate.cpp

@@ -185,7 +185,9 @@ IGL_INLINE bool igl::decimate(
       cost_and_placement(e,V,F,E,EMAP,EF,EI,cost,p);
       C.row(e) = p;
       costs(e) = cost;
-    },10000);
+    },
+    10000
+    );
     for(int e = 0;e<E.rows();e++)
     {
       Q.emplace(costs(e),e,0);

include/igl/dual_contouring.cpp

@@ -0,0 +1,559 @@
+#include "dual_contouring.h"
+#include "quadprog.h"
+#include "parallel_for.h"
+#include <thread>
+#include <mutex>
+#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:
+      using RowVector3S = Eigen::Matrix<Scalar,1,3>;
+      using RowVector4S = Eigen::Matrix<Scalar,1,4>;
+      using Matrix4S = Eigen::Matrix<Scalar,4,4>;
+      using Matrix3S = Eigen::Matrix<Scalar,3,3>;
+      using Vector3S = Eigen::Matrix<Scalar,3,1>;
+      using KeyTriplet = std::tuple<int,int,int>;
+    // 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,[&](const int x)
+        {
+          for(int y = 0;y<ny;y++)
+          {
+            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);
+              }
+            }
+          }
+        },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,
+        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 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,
+  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<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

@@ -0,0 +1,103 @@
+#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);
+  // 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
+  //   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

include/igl/fast_winding_number.cpp

@@ -464,6 +464,8 @@ IGL_INLINE typename Derivedp::Scalar igl::fast_winding_number(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::fast_winding_number<Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(igl::FastWindingNumberBVH const&, float, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 template Eigen::Matrix<float, -1, -1, 0, -1, -1>::Scalar igl::fast_winding_number<Eigen::Matrix<float, -1, -1, 0, -1, -1> >(igl::FastWindingNumberBVH const&, float, Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&);
 // generated by autoexplicit.sh
 template void igl::fast_winding_number<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<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);

include/igl/fit_rotations.cpp

@@ -222,4 +222,5 @@ template void igl::fit_rotations<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen
 template void igl::fit_rotations_planar<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::fit_rotations_planar<Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> >&);
 template void igl::fit_rotations<Eigen::Matrix<float,-1,-1,0,-1,-1>,Eigen::Matrix<float,-1,-1,0,-1,-1> >(Eigen::PlainObjectBase<Eigen::Matrix<float,-1,-1,0,-1,-1> > const &,bool,Eigen::PlainObjectBase<Eigen::Matrix<float,-1,-1,0,-1,-1> > &);
+template void igl::fit_rotations_planar<Eigen::Matrix<double, 2, 2, 0, 2, 2>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 2, 2, 0, 2, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/grid.cpp

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

include/igl/iterative_closest_point.cpp

@@ -8,12 +8,9 @@
 #include "iterative_closest_point.h"
 #include "AABB.h"
 #include "per_face_normals.h"
-#include "matlab_format.h"
 #include "slice.h"
 #include "random_points_on_mesh.h"
 #include "rigid_alignment.h"
-#include "writeOBJ.h"
-#include "writeDMAT.h"
 #include <cassert>
 #include <iostream>
 

include/igl/march_cube.cpp

@@ -126,6 +126,8 @@ IGL_INLINE void igl::march_cube(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+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<int64_t, int, std::hash<int64_t>, std::equal_to<int64_t>, std::allocator<std::pair<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<int64_t, int, std::hash<int64_t>, std::equal_to<int64_t>, std::allocator<std::pair<int64_t const, int> > >&);
 template void igl::march_cube<Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >, double, long, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > 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, -1, 0, -1, -1> >&, long&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, long&, std::unordered_map<int64_t, int, std::hash<int64_t>, std::equal_to<int64_t>, std::allocator<std::pair<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, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(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, -1, 0, -1, -1> >&, unsigned int&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, unsigned int&, std::unordered_map<int64_t, int, std::hash<int64_t>, std::equal_to<int64_t>, std::allocator<std::pair<int64_t const, int> > >&);
 #ifdef WIN32

include/igl/marching_cubes.cpp

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

include/igl/matlab_format.cpp

@@ -34,6 +34,15 @@ IGL_INLINE const Eigen::WithFormat< DerivedM > igl::matlab_format(
     "\n];"));
 }
 
+template <typename DerivedM>
+IGL_INLINE std::string igl::matlab_format_index(
+  const Eigen::MatrixBase<DerivedM> & M,
+  const std::string name)
+{
+  // 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(
@@ -114,6 +123,126 @@ IGL_INLINE const std::string igl::matlab_format(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, -1, 1, 0, 4, 1> > const igl::matlab_format<Eigen::Matrix<double, -1, 1, 0, 4, 1> >(Eigen::DenseBase<Eigen::Matrix<double, -1, 1, 0, 4, 1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, -1, 1, 0, 10, 1> > const igl::matlab_format<Eigen::Matrix<double, -1, 1, 0, 10, 1> >(Eigen::DenseBase<Eigen::Matrix<double, -1, 1, 0, 10, 1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Array<bool, 8, 1, 0, 8, 1> > const igl::matlab_format<Eigen::Array<bool, 8, 1, 0, 8, 1> >(Eigen::DenseBase<Eigen::Array<bool, 8, 1, 0, 8, 1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 2, 8, 0, 2, 8> > const igl::matlab_format<Eigen::Matrix<double, 2, 8, 0, 2, 8> >(Eigen::DenseBase<Eigen::Matrix<double, 2, 8, 0, 2, 8> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 2, 6, 0, 2, 6> > const igl::matlab_format<Eigen::Matrix<double, 2, 6, 0, 2, 6> >(Eigen::DenseBase<Eigen::Matrix<double, 2, 6, 0, 2, 6> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 0, 3, 0, 0, 3> > const igl::matlab_format<Eigen::Matrix<double, 0, 3, 0, 0, 3> >(Eigen::DenseBase<Eigen::Matrix<double, 0, 3, 0, 0, 3> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 0, 2, 0, 0, 2> > const igl::matlab_format<Eigen::Matrix<double, 0, 2, 0, 0, 2> >(Eigen::DenseBase<Eigen::Matrix<double, 0, 2, 0, 0, 2> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 1, 1, 0, 1, 1> > const igl::matlab_format<Eigen::Matrix<double, 1, 1, 0, 1, 1> >(Eigen::DenseBase<Eigen::Matrix<double, 1, 1, 0, 1, 1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Array<bool, -1, 1, 0, -1, 1> > const igl::matlab_format<Eigen::Array<bool, -1, 1, 0, -1, 1> >(Eigen::DenseBase<Eigen::Array<bool, -1, 1, 0, -1, 1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Array<bool, 3, 1, 0, 3, 1> > const igl::matlab_format<Eigen::Array<bool, 3, 1, 0, 3, 1> >(Eigen::DenseBase<Eigen::Array<bool, 3, 1, 0, 3, 1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Array<bool, 2, 1, 0, 2, 1> > const igl::matlab_format<Eigen::Array<bool, 2, 1, 0, 2, 1> >(Eigen::DenseBase<Eigen::Array<bool, 2, 1, 0, 2, 1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 16, 1, 0, 16, 1> > const igl::matlab_format<Eigen::Matrix<double, 16, 1, 0, 16, 1> >(Eigen::DenseBase<Eigen::Matrix<double, 16, 1, 0, 16, 1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 16, 16, 0, 16, 16> > const igl::matlab_format<Eigen::Matrix<double, 16, 16, 0, 16, 16> >(Eigen::DenseBase<Eigen::Matrix<double, 16, 16, 0, 16, 16> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 15, 1, 0, 15, 1> > const igl::matlab_format<Eigen::Matrix<double, 15, 1, 0, 15, 1> >(Eigen::DenseBase<Eigen::Matrix<double, 15, 1, 0, 15, 1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 15, 15, 0, 15, 15> > const igl::matlab_format<Eigen::Matrix<double, 15, 15, 0, 15, 15> >(Eigen::DenseBase<Eigen::Matrix<double, 15, 15, 0, 15, 15> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 14, 1, 0, 14, 1> > const igl::matlab_format<Eigen::Matrix<double, 14, 1, 0, 14, 1> >(Eigen::DenseBase<Eigen::Matrix<double, 14, 1, 0, 14, 1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 14, 14, 0, 14, 14> > const igl::matlab_format<Eigen::Matrix<double, 14, 14, 0, 14, 14> >(Eigen::DenseBase<Eigen::Matrix<double, 14, 14, 0, 14, 14> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 13, 1, 0, 13, 1> > const igl::matlab_format<Eigen::Matrix<double, 13, 1, 0, 13, 1> >(Eigen::DenseBase<Eigen::Matrix<double, 13, 1, 0, 13, 1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 13, 13, 0, 13, 13> > const igl::matlab_format<Eigen::Matrix<double, 13, 13, 0, 13, 13> >(Eigen::DenseBase<Eigen::Matrix<double, 13, 13, 0, 13, 13> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 12, 1, 0, 12, 1> > const igl::matlab_format<Eigen::Matrix<double, 12, 1, 0, 12, 1> >(Eigen::DenseBase<Eigen::Matrix<double, 12, 1, 0, 12, 1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 12, 12, 0, 12, 12> > const igl::matlab_format<Eigen::Matrix<double, 12, 12, 0, 12, 12> >(Eigen::DenseBase<Eigen::Matrix<double, 12, 12, 0, 12, 12> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 11, 1, 0, 11, 1> > const igl::matlab_format<Eigen::Matrix<double, 11, 1, 0, 11, 1> >(Eigen::DenseBase<Eigen::Matrix<double, 11, 1, 0, 11, 1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 11, 11, 0, 11, 11> > const igl::matlab_format<Eigen::Matrix<double, 11, 11, 0, 11, 11> >(Eigen::DenseBase<Eigen::Matrix<double, 11, 11, 0, 11, 11> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 9, 18, 1, 9, 18> > const igl::matlab_format<Eigen::Matrix<double, 9, 18, 1, 9, 18> >(Eigen::DenseBase<Eigen::Matrix<double, 9, 18, 1, 9, 18> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 8, 16, 1, 8, 16> > const igl::matlab_format<Eigen::Matrix<double, 8, 16, 1, 8, 16> >(Eigen::DenseBase<Eigen::Matrix<double, 8, 16, 1, 8, 16> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 80, 160, 1, 80, 160> > const igl::matlab_format<Eigen::Matrix<double, 80, 160, 1, 80, 160> >(Eigen::DenseBase<Eigen::Matrix<double, 80, 160, 1, 80, 160> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 7, 14, 1, 7, 14> > const igl::matlab_format<Eigen::Matrix<double, 7, 14, 1, 7, 14> >(Eigen::DenseBase<Eigen::Matrix<double, 7, 14, 1, 7, 14> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 6, 12, 1, 6, 12> > const igl::matlab_format<Eigen::Matrix<double, 6, 12, 1, 6, 12> >(Eigen::DenseBase<Eigen::Matrix<double, 6, 12, 1, 6, 12> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 5, 10, 1, 5, 10> > const igl::matlab_format<Eigen::Matrix<double, 5, 10, 1, 5, 10> >(Eigen::DenseBase<Eigen::Matrix<double, 5, 10, 1, 5, 10> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 4, 8, 1, 4, 8> > const igl::matlab_format<Eigen::Matrix<double, 4, 8, 1, 4, 8> >(Eigen::DenseBase<Eigen::Matrix<double, 4, 8, 1, 4, 8> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 40, 80, 1, 40, 80> > const igl::matlab_format<Eigen::Matrix<double, 40, 80, 1, 40, 80> >(Eigen::DenseBase<Eigen::Matrix<double, 40, 80, 1, 40, 80> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 3, 6, 1, 3, 6> > const igl::matlab_format<Eigen::Matrix<double, 3, 6, 1, 3, 6> >(Eigen::DenseBase<Eigen::Matrix<double, 3, 6, 1, 3, 6> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 20, 40, 1, 20, 40> > const igl::matlab_format<Eigen::Matrix<double, 20, 40, 1, 20, 40> >(Eigen::DenseBase<Eigen::Matrix<double, 20, 40, 1, 20, 40> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 10, 20, 1, 10, 20> > const igl::matlab_format<Eigen::Matrix<double, 10, 20, 1, 10, 20> >(Eigen::DenseBase<Eigen::Matrix<double, 10, 20, 1, 10, 20> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 9, 9, 0, 9, 9> > const igl::matlab_format<Eigen::Matrix<double, 9, 9, 0, 9, 9> >(Eigen::DenseBase<Eigen::Matrix<double, 9, 9, 0, 9, 9> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 8, 8, 0, 8, 8> > const igl::matlab_format<Eigen::Matrix<double, 8, 8, 0, 8, 8> >(Eigen::DenseBase<Eigen::Matrix<double, 8, 8, 0, 8, 8> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 80, 80, 0, 80, 80> > const igl::matlab_format<Eigen::Matrix<double, 80, 80, 0, 80, 80> >(Eigen::DenseBase<Eigen::Matrix<double, 80, 80, 0, 80, 80> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 80, 1, 0, 80, 1> > const igl::matlab_format<Eigen::Matrix<double, 80, 1, 0, 80, 1> >(Eigen::DenseBase<Eigen::Matrix<double, 80, 1, 0, 80, 1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 7, 7, 0, 7, 7> > const igl::matlab_format<Eigen::Matrix<double, 7, 7, 0, 7, 7> >(Eigen::DenseBase<Eigen::Matrix<double, 7, 7, 0, 7, 7> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 6, 6, 0, 6, 6> > const igl::matlab_format<Eigen::Matrix<double, 6, 6, 0, 6, 6> >(Eigen::DenseBase<Eigen::Matrix<double, 6, 6, 0, 6, 6> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 5, 5, 0, 5, 5> > const igl::matlab_format<Eigen::Matrix<double, 5, 5, 0, 5, 5> >(Eigen::DenseBase<Eigen::Matrix<double, 5, 5, 0, 5, 5> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 40, 40, 0, 40, 40> > const igl::matlab_format<Eigen::Matrix<double, 40, 40, 0, 40, 40> >(Eigen::DenseBase<Eigen::Matrix<double, 40, 40, 0, 40, 40> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 40, 1, 0, 40, 1> > const igl::matlab_format<Eigen::Matrix<double, 40, 1, 0, 40, 1> >(Eigen::DenseBase<Eigen::Matrix<double, 40, 1, 0, 40, 1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 20, 20, 0, 20, 20> > const igl::matlab_format<Eigen::Matrix<double, 20, 20, 0, 20, 20> >(Eigen::DenseBase<Eigen::Matrix<double, 20, 20, 0, 20, 20> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 20, 1, 0, 20, 1> > const igl::matlab_format<Eigen::Matrix<double, 20, 1, 0, 20, 1> >(Eigen::DenseBase<Eigen::Matrix<double, 20, 1, 0, 20, 1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 10, 10, 0, 10, 10> > const igl::matlab_format<Eigen::Matrix<double, 10, 10, 0, 10, 10> >(Eigen::DenseBase<Eigen::Matrix<double, 10, 10, 0, 10, 10> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 2, 4, 1, 2, 4> > const igl::matlab_format<Eigen::Matrix<double, 2, 4, 1, 2, 4> >(Eigen::DenseBase<Eigen::Matrix<double, 2, 4, 1, 2, 4> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 2, 0, 1, 2, 0> > const igl::matlab_format<Eigen::Matrix<double, 2, 0, 1, 2, 0> >(Eigen::DenseBase<Eigen::Matrix<double, 2, 0, 1, 2, 0> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 0, 1, 0, 0, 1> > const igl::matlab_format<Eigen::Matrix<double, 0, 1, 0, 0, 1> >(Eigen::DenseBase<Eigen::Matrix<double, 0, 1, 0, 0, 1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+template Eigen::WithFormat<Eigen::Matrix<double, 2, 0, 0, 2, 0> > const igl::matlab_format<Eigen::Matrix<double, 2, 0, 0, 2, 0> >(Eigen::DenseBase<Eigen::Matrix<double, 2, 0, 0, 2, 0> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+template Eigen::WithFormat<Eigen::Matrix<double, 2, 4, 0, 2, 4> > const igl::matlab_format<Eigen::Matrix<double, 2, 4, 0, 2, 4> >(Eigen::DenseBase<Eigen::Matrix<double, 2, 4, 0, 2, 4> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 9, 1, 0, 9, 1> > const igl::matlab_format<Eigen::Matrix<double, 9, 1, 0, 9, 1> >(Eigen::DenseBase<Eigen::Matrix<double, 9, 1, 0, 9, 1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 8, 1, 0, 8, 1> > const igl::matlab_format<Eigen::Matrix<double, 8, 1, 0, 8, 1> >(Eigen::DenseBase<Eigen::Matrix<double, 8, 1, 0, 8, 1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 7, 1, 0, 7, 1> > const igl::matlab_format<Eigen::Matrix<double, 7, 1, 0, 7, 1> >(Eigen::DenseBase<Eigen::Matrix<double, 7, 1, 0, 7, 1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 6, 1, 0, 6, 1> > const igl::matlab_format<Eigen::Matrix<double, 6, 1, 0, 6, 1> >(Eigen::DenseBase<Eigen::Matrix<double, 6, 1, 0, 6, 1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 5, 1, 0, 5, 1> > const igl::matlab_format<Eigen::Matrix<double, 5, 1, 0, 5, 1> >(Eigen::DenseBase<Eigen::Matrix<double, 5, 1, 0, 5, 1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 4, 1, 0, 4, 1> > const igl::matlab_format<Eigen::Matrix<double, 4, 1, 0, 4, 1> >(Eigen::DenseBase<Eigen::Matrix<double, 4, 1, 0, 4, 1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 10, 1, 0, 10, 1> > const igl::matlab_format<Eigen::Matrix<double, 10, 1, 0, 10, 1> >(Eigen::DenseBase<Eigen::Matrix<double, 10, 1, 0, 10, 1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 99, 99, 0, 99, 99> > const igl::matlab_format<Eigen::Matrix<double, 99, 99, 0, 99, 99> >(Eigen::DenseBase<Eigen::Matrix<double, 99, 99, 0, 99, 99> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::Matrix<double, 99, 1, 0, 99, 1> > const igl::matlab_format<Eigen::Matrix<double, 99, 1, 0, 99, 1> >(Eigen::DenseBase<Eigen::Matrix<double, 99, 1, 0, 99, 1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+template std::basic_string<char, std::char_traits<char>, std::allocator<char> > igl::matlab_format_index<Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+template std::basic_string<char, std::char_traits<char>, std::allocator<char> > igl::matlab_format_index<Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+///////////////////////////////////////////////////
 template Eigen::WithFormat<Eigen::Array<int, -1, 3, 0, -1, 3> > const igl::matlab_format<Eigen::Array<int, -1, 3, 0, -1, 3> >(Eigen::DenseBase<Eigen::Array<int, -1, 3, 0, -1, 3> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
 template Eigen::WithFormat<Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> > const igl::matlab_format<Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> >(Eigen::DenseBase<Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
 template Eigen::WithFormat<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const igl::matlab_format<Eigen::Matrix<int, -1, 3, 1, -1, 3> >(Eigen::DenseBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);

include/igl/matlab_format.h

@@ -44,6 +44,10 @@ namespace igl
   IGL_INLINE const Eigen::WithFormat< DerivedM > matlab_format(
     const Eigen::DenseBase<DerivedM> & M,
     const std::string name = "");
+  template <typename DerivedM>
+  IGL_INLINE std::string matlab_format_index(
+    const Eigen::MatrixBase<DerivedM> & M,
+    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
   // Example:

include/igl/min_quad_with_fixed.cpp

@@ -21,7 +21,8 @@
 #include <unsupported/Eigen/SparseExtra>
 #include <cassert>
 #include <cstdio>
-#include <iostream>
+#include <igl/matlab_format.h>
+#include <type_traits>
 
 template <typename T, typename Derivedknown>
 IGL_INLINE bool igl::min_quad_with_fixed_precompute(
@@ -582,9 +583,312 @@ IGL_INLINE bool igl::min_quad_with_fixed(
   return min_quad_with_fixed_solve(data,B,Y,Beq,Z);
 }
 
+
+template <typename Scalar, int n, int m, bool Hpd>
+IGL_INLINE Eigen::Matrix<Scalar,n,1> igl::min_quad_with_fixed(
+  const Eigen::Matrix<Scalar,n,n> & H,
+  const Eigen::Matrix<Scalar,n,1> & f,
+  const Eigen::Array<bool,n,1> & k,
+  const Eigen::Matrix<Scalar,n,1> & bc,
+  const Eigen::Matrix<Scalar,m,n> & A,
+  const Eigen::Matrix<Scalar,m,1> & b)
+{
+  const auto dyn_n = n == Eigen::Dynamic ? H.rows() : n;
+  const auto dyn_m = m == Eigen::Dynamic ? A.rows() : m;
+  constexpr const int nn = n == Eigen::Dynamic ? Eigen::Dynamic : n+m;
+  const auto dyn_nn = nn == Eigen::Dynamic ? dyn_n+dyn_m : nn;
+  if(dyn_m == 0)
+  {
+    return igl::min_quad_with_fixed<Scalar,n,Hpd>(H,f,k,bc);
+  }
+  // min_x ½ xᵀ H x + xᵀ f   subject to A x = b and x(k) = bc(k)
+  // let zᵀ = [xᵀ λᵀ]
+  // min_z ½ zᵀ [H Aᵀ;A 0] z + zᵀ [f;-b]   z(k) = bc(k)
+  const auto make_HH = [&]()
+  {
+    // Windows can't remember that nn is const.
+    constexpr const int nn = n == Eigen::Dynamic ? Eigen::Dynamic : n+m;
+    Eigen::Matrix<Scalar,nn,nn> HH =
+      Eigen::Matrix<Scalar,nn,nn>::Zero(dyn_nn,dyn_nn);
+    HH.topLeftCorner(dyn_n,dyn_n) = H;
+    HH.bottomLeftCorner(dyn_m,dyn_n) = A;
+    HH.topRightCorner(dyn_n,dyn_m) = A.transpose();
+    return HH;
+  };
+  const Eigen::Matrix<Scalar,nn,nn> HH = make_HH();
+  const auto make_ff  = [&]()
+  {
+    // Windows can't remember that nn is const.
+    constexpr const int nn = n == Eigen::Dynamic ? Eigen::Dynamic : n+m;
+    Eigen::Matrix<Scalar,nn,1> ff(dyn_nn);
+    ff.head(dyn_n) =  f;
+    ff.tail(dyn_m) = -b;
+    return ff;
+  };
+  const Eigen::Matrix<Scalar,nn,1> ff = make_ff();
+  const auto make_kk  = [&]()
+  {
+    // Windows can't remember that nn is const.
+    constexpr const int nn = n == Eigen::Dynamic ? Eigen::Dynamic : n+m;
+    Eigen::Array<bool,nn,1> kk =
+      Eigen::Array<bool,nn,1>::Constant(dyn_nn,1,false);
+    kk.head(dyn_n) =  k;
+    return kk;
+  };
+  const Eigen::Array<bool,nn,1> kk = make_kk();
+  const auto make_bcbc= [&]()
+  {
+    // Windows can't remember that nn is const.
+    constexpr const int nn = n == Eigen::Dynamic ? Eigen::Dynamic : n+m;
+    Eigen::Matrix<Scalar,nn,1> bcbc(dyn_nn);
+    bcbc.head(dyn_n) =  bc;
+    return bcbc;
+  };
+  const Eigen::Matrix<double,nn,1> bcbc = make_bcbc();
+  const Eigen::Matrix<Scalar,nn,1> xx = 
+    min_quad_with_fixed<Scalar,nn,false>(HH,ff,kk,bcbc);
+  return xx.head(dyn_n);
+}
+
+template <typename Scalar, int n, bool Hpd>
+IGL_INLINE Eigen::Matrix<Scalar,n,1> igl::min_quad_with_fixed(
+  const Eigen::Matrix<Scalar,n,n> & H,
+  const Eigen::Matrix<Scalar,n,1> & f,
+  const Eigen::Array<bool,n,1> & k,
+  const Eigen::Matrix<Scalar,n,1> & bc)
+{
+  assert(H.isApprox(H.transpose(),1e-7)); 
+  assert(H.rows() == H.cols());
+  assert(H.rows() == f.size());
+  assert(H.rows() == k.size());
+  assert(H.rows() == bc.size());
+  const auto kcount = k.count();
+  // Everything fixed
+  if(kcount == (Eigen::Dynamic?H.rows():n))
+  {
+    return bc;
+  }
+  // Nothing fixed
+  if(kcount == 0)
+  {
+    // avoid function call
+    typedef Eigen::Matrix<Scalar,n,n> MatrixSn;
+    typedef typename 
+      std::conditional<Hpd,Eigen::LLT<MatrixSn>,Eigen::CompleteOrthogonalDecomposition<MatrixSn>>::type
+      Solver;
+    return Solver(H).solve(-f);
+  }
+  // All-but-one fixed
+  if( (Eigen::Dynamic?H.rows():n)-kcount == 1)
+  {
+    // which one is not fixed?
+    int u = -1;
+    for(int i=0;i<k.size();i++){ if(!k(i)){ u=i; break; } } 
+    assert(u>=0);
+    // min ½ x(u) Huu x(u) + x(u)(fu + H(u,k)bc(k))
+    // Huu x(u) = -(fu + H(u,k) bc(k))
+    // x(u) = (-fu + ∑ -Huj bcj)/Huu
+    Eigen::Matrix<Scalar,n,1> x = bc;
+    x(u) = -f(u);
+    for(int i=0;i<k.size();i++){ if(i!=u){ x(u)-=bc(i)*H(i,u); } } 
+    x(u) /= H(u,u);
+    return x;
+  }
+  // Is there a smart template way to do this?
+  switch(kcount)
+  {
+    case 0: assert(false && "Handled above."); return Eigen::Matrix<Scalar,n,1>();
+    // maxi=16;for i=1:maxi;fprintf('    case %d:\n    {\n     const bool D = (n-%d<=0)||(%d>=n)||(n>%d);\n     return min_quad_with_fixed<Scalar,D?Eigen::Dynamic:n,D?Eigen::Dynamic:%d,Hpd>(H,f,k,bc);\n    }\n',[i i i maxi i]);end
+    case 1:
+    {
+     const bool D = (n-1<=0)||(1>=n)||(n>16);
+     return min_quad_with_fixed<Scalar,D?Eigen::Dynamic:n,D?Eigen::Dynamic:1,Hpd>(H,f,k,bc);
+    }
+    case 2:
+    {
+     const bool D = (n-2<=0)||(2>=n)||(n>16);
+     return min_quad_with_fixed<Scalar,D?Eigen::Dynamic:n,D?Eigen::Dynamic:2,Hpd>(H,f,k,bc);
+    }
+    case 3:
+    {
+     const bool D = (n-3<=0)||(3>=n)||(n>16);
+     return min_quad_with_fixed<Scalar,D?Eigen::Dynamic:n,D?Eigen::Dynamic:3,Hpd>(H,f,k,bc);
+    }
+    case 4:
+    {
+     const bool D = (n-4<=0)||(4>=n)||(n>16);
+     return min_quad_with_fixed<Scalar,D?Eigen::Dynamic:n,D?Eigen::Dynamic:4,Hpd>(H,f,k,bc);
+    }
+    case 5:
+    {
+     const bool D = (n-5<=0)||(5>=n)||(n>16);
+     return min_quad_with_fixed<Scalar,D?Eigen::Dynamic:n,D?Eigen::Dynamic:5,Hpd>(H,f,k,bc);
+    }
+    case 6:
+    {
+     const bool D = (n-6<=0)||(6>=n)||(n>16);
+     return min_quad_with_fixed<Scalar,D?Eigen::Dynamic:n,D?Eigen::Dynamic:6,Hpd>(H,f,k,bc);
+    }
+    case 7:
+    {
+     const bool D = (n-7<=0)||(7>=n)||(n>16);
+     return min_quad_with_fixed<Scalar,D?Eigen::Dynamic:n,D?Eigen::Dynamic:7,Hpd>(H,f,k,bc);
+    }
+    case 8:
+    {
+     const bool D = (n-8<=0)||(8>=n)||(n>16);
+     return min_quad_with_fixed<Scalar,D?Eigen::Dynamic:n,D?Eigen::Dynamic:8,Hpd>(H,f,k,bc);
+    }
+    case 9:
+    {
+     const bool D = (n-9<=0)||(9>=n)||(n>16);
+     return min_quad_with_fixed<Scalar,D?Eigen::Dynamic:n,D?Eigen::Dynamic:9,Hpd>(H,f,k,bc);
+    }
+    case 10:
+    {
+     const bool D = (n-10<=0)||(10>=n)||(n>16);
+     return min_quad_with_fixed<Scalar,D?Eigen::Dynamic:n,D?Eigen::Dynamic:10,Hpd>(H,f,k,bc);
+    }
+    case 11:
+    {
+     const bool D = (n-11<=0)||(11>=n)||(n>16);
+     return min_quad_with_fixed<Scalar,D?Eigen::Dynamic:n,D?Eigen::Dynamic:11,Hpd>(H,f,k,bc);
+    }
+    case 12:
+    {
+     const bool D = (n-12<=0)||(12>=n)||(n>16);
+     return min_quad_with_fixed<Scalar,D?Eigen::Dynamic:n,D?Eigen::Dynamic:12,Hpd>(H,f,k,bc);
+    }
+    case 13:
+    {
+     const bool D = (n-13<=0)||(13>=n)||(n>16);
+     return min_quad_with_fixed<Scalar,D?Eigen::Dynamic:n,D?Eigen::Dynamic:13,Hpd>(H,f,k,bc);
+    }
+    case 14:
+    {
+     const bool D = (n-14<=0)||(14>=n)||(n>16);
+     return min_quad_with_fixed<Scalar,D?Eigen::Dynamic:n,D?Eigen::Dynamic:14,Hpd>(H,f,k,bc);
+    }
+    case 15:
+    {
+     const bool D = (n-15<=0)||(15>=n)||(n>16);
+     return min_quad_with_fixed<Scalar,D?Eigen::Dynamic:n,D?Eigen::Dynamic:15,Hpd>(H,f,k,bc);
+    }
+    case 16:
+    {
+     const bool D = (n-16<=0)||(16>=n)||(n>16);
+     return min_quad_with_fixed<Scalar,D?Eigen::Dynamic:n,D?Eigen::Dynamic:16,Hpd>(H,f,k,bc);
+    }
+    default:
+      return min_quad_with_fixed<Scalar,Eigen::Dynamic,Eigen::Dynamic,Hpd>(H,f,k,bc);
+  }
+}
+
+template <typename Scalar, int n, int kcount, bool Hpd>
+IGL_INLINE Eigen::Matrix<Scalar,n,1> igl::min_quad_with_fixed(
+  const Eigen::Matrix<Scalar,n,n> & H,
+  const Eigen::Matrix<Scalar,n,1> & f,
+  const Eigen::Array<bool,n,1> & k,
+  const Eigen::Matrix<Scalar,n,1> & bc)
+{
+  // 0 and n should be handle outside this function
+  static_assert(kcount==Eigen::Dynamic || kcount>0                  ,"");
+  static_assert(kcount==Eigen::Dynamic || kcount<n                  ,"");
+  const int ucount = n==Eigen::Dynamic ? Eigen::Dynamic : n-kcount;
+  static_assert(kcount==Eigen::Dynamic || ucount+kcount == n        ,"");
+  static_assert((n==Eigen::Dynamic) == (ucount==Eigen::Dynamic),"");
+  static_assert((kcount==Eigen::Dynamic) == (ucount==Eigen::Dynamic),"");
+  assert((n==Eigen::Dynamic) || n == H.rows());
+  assert((kcount==Eigen::Dynamic) || kcount == k.count());
+  typedef Eigen::Matrix<Scalar,ucount,ucount> MatrixSuu;
+  typedef Eigen::Matrix<Scalar,ucount,kcount> MatrixSuk;
+  typedef Eigen::Matrix<Scalar,n,1>      VectorSn;
+  typedef Eigen::Matrix<Scalar,ucount,1> VectorSu;
+  typedef Eigen::Matrix<Scalar,kcount,1> VectorSk;
+  const auto dyn_n = n==Eigen::Dynamic ? H.rows() : n;
+  const auto dyn_kcount = kcount==Eigen::Dynamic ? k.count() : kcount;
+  const auto dyn_ucount = ucount==Eigen::Dynamic ? dyn_n- dyn_kcount : ucount;
+  // For ucount==2 or kcount==2 this calls the coefficient initiliazer rather
+  // than the size initilizer, but I guess that's ok.
+  MatrixSuu Huu(dyn_ucount,dyn_ucount);
+  MatrixSuk Huk(dyn_ucount,dyn_kcount);
+  VectorSu mrhs(dyn_ucount);
+  VectorSk  bck(dyn_kcount);
+  {
+    int ui = 0;
+    int ki = 0;
+    for(int i = 0;i<dyn_n;i++)
+    {
+      if(k(i))
+      {
+        bck(ki) = bc(i);
+        ki++;
+      }else
+      {
+        mrhs(ui) = f(i);
+        int uj = 0;
+        int kj = 0;
+        for(int j = 0;j<dyn_n;j++)
+        {
+          if(k(j))
+          {
+            Huk(ui,kj) = H(i,j);
+            kj++;
+          }else
+          {
+            Huu(ui,uj) = H(i,j);
+            uj++;
+          }
+        }
+        ui++;
+      }
+    }
+  }
+  mrhs += Huk * bck;
+  typedef typename 
+    std::conditional<Hpd,
+      Eigen::LLT<MatrixSuu>,
+      // LDLT should be faster for indefinite problems but already found some
+      // cases where it was too inaccurate when called via quadprog_primal.
+      // Ideally this function takes LLT,LDLT, or
+      // CompleteOrthogonalDecomposition as a template parameter.  "template
+      // template" parameters did work because LLT,LDLT have different number of
+      // template parameters from CompleteOrthogonalDecomposition.  Perhaps
+      // there's a way to take advantage of LLT and LDLT's default template
+      // parameters (I couldn't figure out how).
+      Eigen::CompleteOrthogonalDecomposition<MatrixSuu>>::type
+    Solver;
+  VectorSu xu = Solver(Huu).solve(-mrhs);
+  VectorSn x(dyn_n);
+  {
+    int ui = 0;
+    int ki = 0;
+    for(int i = 0;i<dyn_n;i++)
+    {
+      if(k(i))
+      {
+        x(i) = bck(ki);
+        ki++;
+      }else
+      {
+        x(i) = xu(ui);
+        ui++;
+      }
+    }
+  }
+  return x;
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template Eigen::Matrix<double, -1, 1, ((Eigen::StorageOptions)0) | ((((-1) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((-1) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), -1, 1> igl::min_quad_with_fixed<double, -1, -1, true>(Eigen::Matrix<double, -1, -1, ((Eigen::StorageOptions)0) | ((((-1) == (1)) && ((-1) != (1))) ? ((Eigen::StorageOptions)1) : ((((-1) == (1)) && ((-1) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), -1, -1> const&, Eigen::Matrix<double, -1, 1, ((Eigen::StorageOptions)0) | ((((-1) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((-1) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), -1, 1> const&, Eigen::Array<bool, -1, 1, ((Eigen::StorageOptions)0) | ((((-1) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((-1) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), -1, 1> const&, Eigen::Matrix<double, -1, 1, ((Eigen::StorageOptions)0) | ((((-1) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((-1) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), -1, 1> const&, Eigen::Matrix<double, -1, -1, ((Eigen::StorageOptions)0) | ((((-1) == (1)) && ((-1) != (1))) ? ((Eigen::StorageOptions)1) : ((((-1) == (1)) && ((-1) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), -1, -1> const&, Eigen::Matrix<double, -1, 1, ((Eigen::StorageOptions)0) | ((((-1) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((-1) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), -1, 1> const&);
+// generated by autoexplicit.sh
+template Eigen::Matrix<double, 8, 1, ((Eigen::StorageOptions)0) | ((((8) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((8) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 8, 1> igl::min_quad_with_fixed<double, 8, 2, true>(Eigen::Matrix<double, 8, 8, ((Eigen::StorageOptions)0) | ((((8) == (1)) && ((8) != (1))) ? ((Eigen::StorageOptions)1) : ((((8) == (1)) && ((8) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 8, 8> const&, Eigen::Matrix<double, 8, 1, ((Eigen::StorageOptions)0) | ((((8) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((8) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 8, 1> const&, Eigen::Array<bool, 8, 1, ((Eigen::StorageOptions)0) | ((((8) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((8) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 8, 1> const&, Eigen::Matrix<double, 8, 1, ((Eigen::StorageOptions)0) | ((((8) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((8) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 8, 1> const&, Eigen::Matrix<double, 2, 8, ((Eigen::StorageOptions)0) | ((((2) == (1)) && ((8) != (1))) ? ((Eigen::StorageOptions)1) : ((((8) == (1)) && ((2) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 2, 8> const&, Eigen::Matrix<double, 2, 1, ((Eigen::StorageOptions)0) | ((((2) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((2) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 2, 1> const&);
+// generated by autoexplicit.sh
+template Eigen::Matrix<double, 3, 1, ((Eigen::StorageOptions)0) | ((((3) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 3, 1> igl::min_quad_with_fixed<double, 3, 1, true>(Eigen::Matrix<double, 3, 3, ((Eigen::StorageOptions)0) | ((((3) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)1) : ((((3) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 3, 3> const&, Eigen::Matrix<double, 3, 1, ((Eigen::StorageOptions)0) | ((((3) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 3, 1> const&, Eigen::Array<bool, 3, 1, ((Eigen::StorageOptions)0) | ((((3) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 3, 1> const&, Eigen::Matrix<double, 3, 1, ((Eigen::StorageOptions)0) | ((((3) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 3, 1> const&, Eigen::Matrix<double, 1, 3, ((Eigen::StorageOptions)0) | ((((1) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)1) : ((((3) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 1, 3> const&, Eigen::Matrix<double, 1, 1, ((Eigen::StorageOptions)0) | ((((1) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 1, 1> const&);
+template Eigen::Matrix<double, 2, 1, ((Eigen::StorageOptions)0) | ((((2) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((2) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 2, 1> igl::min_quad_with_fixed<double, 2, 0, true>(Eigen::Matrix<double, 2, 2, ((Eigen::StorageOptions)0) | ((((2) == (1)) && ((2) != (1))) ? ((Eigen::StorageOptions)1) : ((((2) == (1)) && ((2) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 2, 2> const&, Eigen::Matrix<double, 2, 1, ((Eigen::StorageOptions)0) | ((((2) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((2) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 2, 1> const&, Eigen::Array<bool, 2, 1, ((Eigen::StorageOptions)0) | ((((2) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((2) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 2, 1> const&, Eigen::Matrix<double, 2, 1, ((Eigen::StorageOptions)0) | ((((2) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((2) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 2, 1> const&, Eigen::Matrix<double, 0, 2, ((Eigen::StorageOptions)0) | ((((0) == (1)) && ((2) != (1))) ? ((Eigen::StorageOptions)1) : ((((2) == (1)) && ((0) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 0, 2> const&, Eigen::Matrix<double, 0, 1, ((Eigen::StorageOptions)0) | ((((0) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((0) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 0, 1> const&);
+template Eigen::Matrix<double, 3, 1, ((Eigen::StorageOptions)0) | ((((3) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 3, 1> igl::min_quad_with_fixed<double, 3, 0, true>(Eigen::Matrix<double, 3, 3, ((Eigen::StorageOptions)0) | ((((3) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)1) : ((((3) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 3, 3> const&, Eigen::Matrix<double, 3, 1, ((Eigen::StorageOptions)0) | ((((3) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 3, 1> const&, Eigen::Array<bool, 3, 1, ((Eigen::StorageOptions)0) | ((((3) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 3, 1> const&, Eigen::Matrix<double, 3, 1, ((Eigen::StorageOptions)0) | ((((3) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 3, 1> const&, Eigen::Matrix<double, 0, 3, ((Eigen::StorageOptions)0) | ((((0) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)1) : ((((3) == (1)) && ((0) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 0, 3> const&, Eigen::Matrix<double, 0, 1, ((Eigen::StorageOptions)0) | ((((0) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((0) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 0, 1> const&);
 template bool igl::min_quad_with_fixed<double, Eigen::Matrix<int, -1, 1, 0, -1, 1>, 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<double, -1, 1, 0, -1, 1> >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 #if EIGEN_VERSION_AT_LEAST(3,3,0)
 #else

include/igl/min_quad_with_fixed.h

@@ -117,6 +117,52 @@ namespace igl
     const Eigen::MatrixBase<DerivedBeq> & Beq,
     const bool pd,
     Eigen::PlainObjectBase<DerivedZ> & Z);
+
+  // Dense version optimized for very small, known at compile time sizes. Still
+  // works for Eigen::Dynamic (and then everything needs to be Dynamic).
+  //
+  // min_x ½ xᵀ H x + xᵀ f  
+  // subject to 
+  //   A x = b
+  //   x(i) = bc(i) iff k(i)==true
+  //
+  // Templates:
+  //   Scalar  (e.g., double)
+  //   n  #H or Eigen::Dynamic if not known at compile time
+  //   m  #A or Eigen::Dynamic if not known at compile time
+  //   Hpd  whether H is positive definite (LLT used) or not (QR used)
+  // Inputs:
+  //   H  #H by #H quadratic coefficients (only lower triangle used)
+  //   f  #H linear coefficients
+  //   k  #H list of flags whether to fix value
+  //   bc  #H value to fix to (if !k(i) then bc(i) is ignored)
+  //   A  #A by #H list of linear equality constraint coefficients, must be
+  //     linearly independent (with self and fixed value constraints)
+  //   b  #A list of linear equality right-hand sides
+  // Returns #H-long solution x
+  template <typename Scalar, int n, int m, bool Hpd=true>
+  IGL_INLINE Eigen::Matrix<Scalar,n,1> min_quad_with_fixed(
+    const Eigen::Matrix<Scalar,n,n> & H,
+    const Eigen::Matrix<Scalar,n,1> & f,
+    const Eigen::Array<bool,n,1> & k,
+    const Eigen::Matrix<Scalar,n,1> & bc,
+    const Eigen::Matrix<Scalar,m,n> & A,
+    const Eigen::Matrix<Scalar,m,1> & b);
+  template <typename Scalar, int n, bool Hpd=true>
+  IGL_INLINE Eigen::Matrix<Scalar,n,1> min_quad_with_fixed(
+    const Eigen::Matrix<Scalar,n,n> & H,
+    const Eigen::Matrix<Scalar,n,1> & f,
+    const Eigen::Array<bool,n,1> & k,
+    const Eigen::Matrix<Scalar,n,1> & bc);
+  // Special wrapper where the number of constrained values (i.e., true values
+  // in k) is exposed as a template parameter. Not intended to be called
+  // directly. The overhead of calling the overloads above is already minimal.
+  template <typename Scalar, int n, int kcount, bool Hpd/*no default*/>
+  IGL_INLINE Eigen::Matrix<Scalar,n,1> min_quad_with_fixed(
+    const Eigen::Matrix<Scalar,n,n> & H,
+    const Eigen::Matrix<Scalar,n,1> & f,
+    const Eigen::Array<bool,n,1> & k,
+    const Eigen::Matrix<Scalar,n,1> & bc);
 }
 
 template <typename T>

include/igl/offset_surface.cpp

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

include/igl/per_corner_normals.cpp

@@ -10,6 +10,7 @@
 #include "vertex_triangle_adjacency.h"
 #include "per_face_normals.h"
 #include "PI.h"
+#include "doublearea.h"
 
 template <typename DerivedV, typename DerivedF, typename DerivedCN>
 IGL_INLINE void igl::per_corner_normals(
@@ -101,9 +102,101 @@ IGL_INLINE void igl::per_corner_normals(
     }
   }
 }
+
+
+template <
+  typename DerivedV, 
+  typename DerivedI, 
+  typename DerivedC, 
+  typename DerivedN,
+  typename DerivedVV,
+  typename DerivedFF,
+  typename DerivedJ,
+  typename DerivedNN>
+IGL_INLINE void igl::per_corner_normals(
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedI> & I,
+  const Eigen::MatrixBase<DerivedC> & C,
+  const typename DerivedV::Scalar corner_threshold,
+  Eigen::PlainObjectBase<DerivedN>  & N,
+  Eigen::PlainObjectBase<DerivedVV> & VV,
+  Eigen::PlainObjectBase<DerivedFF> & FF,
+  Eigen::PlainObjectBase<DerivedJ>  & J,
+  Eigen::PlainObjectBase<DerivedNN> & NN)
+{
+  const Eigen::Index m = C.size()-1;
+  typedef Eigen::Index Index;
+  Eigen::MatrixXd FN;
+  per_face_normals(V,I,C,FN,VV,FF,J);
+  typedef typename DerivedN::Scalar Scalar;
+  Eigen::Matrix<Scalar,Eigen::Dynamic,1> AA;
+  doublearea(VV,FF,AA);
+  // VF[i](j) = p means p is the jth face incident on vertex i
+  // to-do micro-optimization to avoid vector<vector>
+  std::vector<std::vector<Eigen::Index>> VF(V.rows());
+  for(Eigen::Index p = 0;p<m;p++)
+  {
+    // number of faces/vertices in this simple polygon
+    const Index np = C(p+1)-C(p);
+    for(Eigen::Index i = 0;i<np;i++)
+    {
+      VF[I(C(p)+i)].push_back(p);
+    }
+  }
+
+  N.resize(I.rows(),3);
+  for(Eigen::Index p = 0;p<m;p++)
+  {
+    // number of faces/vertices in this simple polygon
+    const Index np = C(p+1)-C(p);
+    Eigen::Matrix<Scalar,3,1> fn = FN.row(p);
+    for(Eigen::Index i = 0;i<np;i++)
+    {
+      N.row(C(p)+i).setZero();
+      // Loop over faces sharing this vertex
+      for(const auto & n : VF[I(C(p)+i)])
+      {
+        Eigen::Matrix<Scalar,3,1> ifn = FN.row(n);
+        // dot product between face's normal and other face's normal
+        Scalar dp = fn.dot(ifn);
+        if(dp > cos(corner_threshold*igl::PI/180))
+        {
+          // add to running sum
+          N.row(C(p)+i) += AA(n) * ifn;
+        }
+      }
+      N.row(C(p)+i).normalize();
+    }
+  }
+
+  // Relies on order of FF 
+  NN.resize(FF.rows()*3,3);
+  {
+    Eigen::Index k = 0;
+    for(Eigen::Index p = 0;p<m;p++)
+    {
+      // number of faces/vertices in this simple polygon
+      const Index np = C(p+1)-C(p);
+      for(Eigen::Index i = 0;i<np;i++)
+      {
+        assert(FF(k,0) == I(C(p)+((i+0)%np)));
+        assert(FF(k,1) == I(C(p)+((i+1)%np)));
+        assert(FF(k,2) == V.rows()+p);
+        NN.row(k*3+0) = N.row(C(p)+((i+0)%np));
+        NN.row(k*3+1) = N.row(C(p)+((i+1)%np));
+        NN.row(k*3+2) = FN.row(p);
+        k++;
+      }
+    }
+    assert(k == FF.rows());
+  }
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::per_corner_normals<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -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>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, 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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
 template void igl::per_corner_normals<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, double, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::per_corner_normals<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<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, double, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::per_corner_normals<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::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, double, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);

include/igl/per_corner_normals.h

@@ -55,6 +55,38 @@ namespace igl
     const std::vector<std::vector<IndexType> >& VF,
     const double corner_threshold,
     Eigen::PlainObjectBase<DerivedCN> & CN);
+  // Inputs:
+  //   V  #V by 3 list of mesh vertex positions
+  //   I  #I vectorized list of polygon corner indices into rows of some matrix V
+  //   C  #polygons+1 list of cumulative polygon sizes so that C(i+1)-C(i) = size of
+  //     the ith polygon, and so I(C(i)) through I(C(i+1)-1) are the indices of
+  //     the ith polygon
+  //   corner_threshold  threshold in degrees on sharp angles
+  // Outputs:
+  //   N  #I by 3 list of per corner normals
+  //   VV  #I+#polygons by 3 list of auxiliary triangle mesh vertex positions
+  //   FF  #I by 3 list of triangle indices into rows of VV
+  //   J  #I list of indices into original polygons
+  //   NN  #FF by 3 list of normals for each auxiliary triangle
+  template <
+    typename DerivedV, 
+    typename DerivedI, 
+    typename DerivedC, 
+    typename DerivedN,
+    typename DerivedVV,
+    typename DerivedFF,
+    typename DerivedJ,
+    typename DerivedNN>
+  IGL_INLINE void per_corner_normals(
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedI> & I,
+    const Eigen::MatrixBase<DerivedC> & C,
+    const typename DerivedV::Scalar corner_threshold,
+    Eigen::PlainObjectBase<DerivedN>  & N,
+    Eigen::PlainObjectBase<DerivedVV> & VV,
+    Eigen::PlainObjectBase<DerivedFF> & FF,
+    Eigen::PlainObjectBase<DerivedJ>  & J,
+    Eigen::PlainObjectBase<DerivedNN> & NN);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/per_face_normals.cpp

@@ -103,21 +103,78 @@ IGL_INLINE void igl::per_face_normals_stable(
 
 }
 
+#include "cotmatrix.h"
+
+template <
+  typename DerivedV,
+  typename DerivedI,
+  typename DerivedC,
+  typename DerivedN,
+  typename DerivedVV,
+  typename DerivedFF,
+  typename DerivedJ>
+IGL_INLINE void igl::per_face_normals(
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedI> & I,
+  const Eigen::MatrixBase<DerivedC> & C,
+  Eigen::PlainObjectBase<DerivedN> & N,
+  Eigen::PlainObjectBase<DerivedVV> & VV,
+  Eigen::PlainObjectBase<DerivedFF> & FF,
+  Eigen::PlainObjectBase<DerivedJ> & J)
+{
+  assert(V.cols() == 3);
+  typedef Eigen::Index Index;
+  typedef typename DerivedN::Scalar Scalar;
+  // Use Bunge et al. algorithm in igl::cotmatrix to insert a point for each
+  // polygon which minimizes squared area.
+  {
+    Eigen::SparseMatrix<Scalar> _1,_2,P;
+    igl::cotmatrix(V,I,C,_1,_2,P);
+    VV = P*V;
+  }
+  // number of polygons
+  const Eigen::Index m = C.size()-1;
+  N.resize(m,3);
+  FF.resize(C(m),3);
+  J.resize(C(m));
+  {
+    Eigen::Index k = 0;
+    for(Eigen::Index p = 0;p<m;p++)
+    {
+      N.row(p).setZero();
+      // number of faces/vertices in this simple polygon
+      const Index np = C(p+1)-C(p);
+      for(Eigen::Index i = 0;i<np;i++)
+      {
+        FF.row(k) << 
+          I(C(p)+((i+0)%np)),
+          I(C(p)+((i+1)%np)),
+          V.rows()+p;
+        J(k) = p;
+        k++;
+        typedef Eigen::Matrix<Scalar,1,3> V3;
+        N.row(p) +=
+          V3(VV.row(I(C(p)+((i+0)%np)))-VV.row(V.rows()+p)).cross(
+          V3(VV.row(I(C(p)+((i+1)%np)))-VV.row(V.rows()+p)));
+      }
+      // normalize to take average
+      N.row(p) /= N.row(p).stableNorm();
+    }
+    assert(k == FF.rows());
+  }
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::per_face_normals<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -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>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, 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<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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::per_face_normals<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<float, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> >&);
-// generated by autoexplicit.sh
 template void igl::per_face_normals<Eigen::Matrix<float, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<float, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> >&);
-// generated by autoexplicit.sh
 template void igl::per_face_normals<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::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, -1, 0, -1, -1> >&);
-// Nonsense template
+// Nonsense template. Where'd this come from? AABB nonsense?
 namespace igl{template<> void per_face_normals<Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&){} }
-// generated by autoexplicit.sh
 template void igl::per_face_normals<Eigen::Matrix<float, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<float, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> >&);
-// generated by autoexplicit.sh
 template void igl::per_face_normals<Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> >&);
-// generated by autoexplicit.sh
 template void igl::per_face_normals<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<float, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> >&);
 template void igl::per_face_normals<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > 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::per_face_normals<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<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);

include/igl/per_face_normals.h

@@ -42,6 +42,34 @@ namespace igl
     const Eigen::MatrixBase<DerivedV>& V,
     const Eigen::MatrixBase<DerivedF>& F,
     Eigen::PlainObjectBase<DerivedN> & N);
+  // Inputs:
+  //   V  #V by 3 list of mesh vertex positions
+  //   I  #I vectorized list of polygon corner indices into rows of some matrix V
+  //   C  #polygons+1 list of cumulative polygon sizes so that C(i+1)-C(i) = size of
+  //     the ith polygon, and so I(C(i)) through I(C(i+1)-1) are the indices of
+  //     the ith polygon
+  //   corner_threshold  threshold in degrees on sharp angles
+  // Outputs:
+  //   N  #polygons by 3 list of per face normals
+  //   VV  #I+#polygons by 3 list of auxiliary triangle mesh vertex positions
+  //   FF  #I by 3 list of triangle indices into rows of VV
+  //   J  #I list of indices into original polygons
+  template <
+    typename DerivedV,
+    typename DerivedI,
+    typename DerivedC,
+    typename DerivedN,
+    typename DerivedVV,
+    typename DerivedFF,
+    typename DerivedJ>
+  IGL_INLINE void per_face_normals(
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedI> & I,
+    const Eigen::MatrixBase<DerivedC> & C,
+    Eigen::PlainObjectBase<DerivedN> & N,
+    Eigen::PlainObjectBase<DerivedVV> & VV,
+    Eigen::PlainObjectBase<DerivedFF> & FF,
+    Eigen::PlainObjectBase<DerivedJ> & J);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/polygon_corners.cpp

@@ -34,8 +34,37 @@ IGL_INLINE void igl::polygon_corners(
   I = Eigen::Map<DerivedI>(vI.data(),vI.size());
 }
 
+template <
+  typename DerivedQ, 
+  typename DerivedI,
+  typename DerivedC>
+IGL_INLINE void igl::polygon_corners(
+  const Eigen::MatrixBase<DerivedQ> & Q,
+  Eigen::PlainObjectBase<DerivedI> & I,
+  Eigen::PlainObjectBase<DerivedC> & C)
+{
+  I.resize(Q.size());
+  C.resize(Q.rows()+1);
+  Eigen::Index c = 0;
+  C(0) = 0;
+  for(Eigen::Index p = 0;p<Q.rows();p++)
+  {
+    Eigen::Index np = 0;
+    for(Eigen::Index i = 0;i<Q.cols();i++)
+    {
+      if(Q(p,i) == -1){ break;}
+      I(c++) = Q(p,i);
+      np++;
+    }
+    C(p+1) = C(p)+np;
+  }
+  I.conservativeResize(c);
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::polygon_corners<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::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+// generated by autoexplicit.sh
 template void igl::polygon_corners<int, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/polygon_corners.h

@@ -32,6 +32,20 @@ namespace igl
     const std::vector<std::vector<PType> > & P,
     Eigen::PlainObjectBase<DerivedI> & I,
     Eigen::PlainObjectBase<DerivedC> & C);
+  // Convert a pure k-gon list of polygon mesh indices to list of polygon corners
+  // and sizes
+  //
+  // Inputs:
+  //   Q  #Q by k list of polygon indices (ith row is a k-gon, unless Q(i,j) =
+  //     -1 then it's a j-gon)
+  template <
+    typename DerivedQ, 
+    typename DerivedI,
+    typename DerivedC>
+  IGL_INLINE void polygon_corners(
+    const Eigen::MatrixBase<DerivedQ> & Q,
+    Eigen::PlainObjectBase<DerivedI> & I,
+    Eigen::PlainObjectBase<DerivedC> & C);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/quadprog.cpp

@@ -0,0 +1,124 @@
+#include "quadprog.h"
+#include "min_quad_with_fixed.h"
+#include <igl/matlab_format.h>
+#include <iostream>
+
+template <typename Scalar, int n, int m>
+IGL_INLINE Eigen::Matrix<Scalar,n,1> igl::quadprog(
+  const Eigen::Matrix<Scalar,n,n> & H,
+  const Eigen::Matrix<Scalar,n,1> & f,
+  const Eigen::Matrix<Scalar,m,n> & A,
+  const Eigen::Matrix<Scalar,m,1> & b,
+  const Eigen::Matrix<Scalar,n,1> & lb,
+  const Eigen::Matrix<Scalar,n,1> & ub)
+{
+  // Alec 16/2/2021:
+  // igl::quadprog implements a very simple primal active set method. The new
+  // igl::min_quad_with_fixed is very fast for small dense problems so the
+  // iterations of igl::quadprog become very fast. Even if it ends up doing many
+  // more iterations than igl::copyleft::quadprog it would be much faster (in
+  // reality it doesn't do that many more iterations). It's a healthy 10-100x
+  // faster than igl::copyleft::quadprog for specific cases of QPs.
+  //
+  // Unfortunately, that set is limited. igl::quadprog is really good at tiny
+  // box-constrained QPs with a positive definite objective (like the kind that show
+  // up in dual contouring). igl::copyleft::quadprog handles more general problems
+  // (and also starts to beat igl::quadprog when the number of variables gets over
+  // ~20). I tried extending igl::quadprog so that we could use it for
+  // igl::copyleft::progressive_hulls and drop igl::copyleft::quadprog but it was
+  // trickier than I thought. Something like qpmad or the non GPL version of
+  // quadrog++ would be good future PR.
+  //
+  typedef Eigen::Matrix<Scalar,n,1> VectorSn;
+  typedef Eigen::Array<bool,n,1>    Arraybn;
+  assert( (lb.array() < ub.array() ).all() );
+  const int dyn_n = n == Eigen::Dynamic ? H.rows() : n;
+  VectorSn x(dyn_n);
+  VectorSn bc = VectorSn::Constant(dyn_n,1,-1e26);
+  Arraybn k = Arraybn::Constant(dyn_n,1,false);
+  Eigen::Index iter;
+  // n³ is probably way too conservative. 
+  for(iter = 0;iter<dyn_n*dyn_n*dyn_n;iter++)
+  {
+    // For dual contouring 99% of the time the active set is empty.
+    // Optimize for this common case.
+    // Windows needs template arguments spelled out
+    x = min_quad_with_fixed<Scalar,n,m>(H,f,k,bc,A,b);
+    // constraint violations 
+    VectorSn vl = lb-x;
+    VectorSn vu = x-ub;
+
+    // try to add/remove constraints
+    Eigen::Index best_add = -1; Scalar worst_offense = 0;
+    bool add_lower;
+    Eigen::Index best_remove = -1; Scalar worst_lambda = 0;
+    for(Eigen::Index i = 0;i<dyn_n;i++)
+    {
+      if(vl(i)>worst_offense)
+      {
+        best_add = i;
+        add_lower = true;
+        worst_offense = vl(i);
+      }
+      if(vu(i)>worst_offense)
+      {
+        best_add = i;
+        add_lower = false;
+        worst_offense = vu(i);
+      }
+      // bias toward adding constraints
+      if(best_add<0 && k(i))
+      {
+        const Scalar sign = bc(i)==ub(i)?1:-1;
+        const Scalar lambda_i = sign * (H.row(i)*x+f(i));
+        if(lambda_i > worst_lambda)
+        {
+          best_remove = i;
+          worst_lambda = lambda_i;
+        }
+      }
+    }
+    // bias toward adding constraints
+    if(best_add >= 0)
+    {
+      const auto i = best_add;
+      assert(!k(i));
+      bc(i) = add_lower ? lb(i) : ub(i);
+      k(i) = true;
+    }else if(best_remove >= 0)
+    {
+      const auto i = best_remove;
+      assert(k(i));
+      k(i) = false;
+    }else /*if(best_add < 0 && best_remove < 0)*/
+    {
+      return x;
+    }
+  }
+  // Should never happen.
+  assert(false && "quadprog failed after too many iterations");
+  return VectorSn::Zero(dyn_n);
+}
+
+template <typename Scalar, int n>
+IGL_INLINE Eigen::Matrix<Scalar,n,1> igl::quadprog(
+  const Eigen::Matrix<Scalar,n,n> & H,
+  const Eigen::Matrix<Scalar,n,1> & f,
+  const Eigen::Matrix<Scalar,n,1> & lb,
+  const Eigen::Matrix<Scalar,n,1> & ub)
+{
+  const int m = n == Eigen::Dynamic ? Eigen::Dynamic : 0;
+  // Windows needs template parameters spelled out
+  return quadprog<Scalar,n,m>(
+    H,f,
+    Eigen::Matrix<Scalar,m,n>(0,H.cols()),
+    Eigen::Matrix<Scalar,m,1>(0,1),
+    lb,ub);
+}
+
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template Eigen::Matrix<double, 2, 1, ((Eigen::StorageOptions)0) | ((((2) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((2) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 2, 1> igl::quadprog<double, 2>(Eigen::Matrix<double, 2, 2, ((Eigen::StorageOptions)0) | ((((2) == (1)) && ((2) != (1))) ? ((Eigen::StorageOptions)1) : ((((2) == (1)) && ((2) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 2, 2> const&, Eigen::Matrix<double, 2, 1, ((Eigen::StorageOptions)0) | ((((2) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((2) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 2, 1> const&, Eigen::Matrix<double, 2, 1, ((Eigen::StorageOptions)0) | ((((2) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((2) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 2, 1> const&, Eigen::Matrix<double, 2, 1, ((Eigen::StorageOptions)0) | ((((2) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((2) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 2, 1> const&);
+template Eigen::Matrix<double, 3, 1, ((Eigen::StorageOptions)0) | ((((3) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 3, 1> igl::quadprog<double, 3>(Eigen::Matrix<double, 3, 3, ((Eigen::StorageOptions)0) | ((((3) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)1) : ((((3) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 3, 3> const&, Eigen::Matrix<double, 3, 1, ((Eigen::StorageOptions)0) | ((((3) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 3, 1> const&, Eigen::Matrix<double, 3, 1, ((Eigen::StorageOptions)0) | ((((3) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 3, 1> const&, Eigen::Matrix<double, 3, 1, ((Eigen::StorageOptions)0) | ((((3) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 3, 1> const&);
+#endif

include/igl/quadprog.h

@@ -0,0 +1,87 @@
+#ifndef IGL_QUADPROG_H
+#define IGL_QUADPROG_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Solve a convex quadratic program. Optimized for small dense problems.
+  // Still works for Eigen::Dynamic (and then everything needs to be Dynamic).
+  //
+  // min_x ½ xᵀ H x + xᵀf
+  // subject to:
+  //   lbi ≤ Ai x ≤ ubi
+  //   lb ≤ x ≤ u
+  //
+  // Templates:
+  //   Scalar  (e.g., double)
+  //   n  #H or Eigen::Dynamic if not known at compile time
+  //   ni  #Ai or Eigen::Dynamic if not known at compile time
+  // Inputs:
+  //   H  #H by #H quadratic coefficients (only lower triangle used)
+  //   f  #H linear coefficients
+  //   Ai  #Ai by #H list of linear equality constraint coefficients
+  //   lbi  #Ai list of linear equality lower bounds
+  //   ubi  #Ai list of linear equality upper bounds
+  //   lb  #H list of lower bounds
+  //   ub  #H list of lower bounds
+  // Returns #H-long solution x
+  template <typename Scalar, int n, int ni>
+  IGL_INLINE Eigen::Matrix<Scalar,n,1> quadprog(
+    const Eigen::Matrix<Scalar,n,n> & H,
+    const Eigen::Matrix<Scalar,n,1> & f,
+    const Eigen::Matrix<Scalar,ni,n> & Ai,
+    const Eigen::Matrix<Scalar,ni,1> & lbi,
+    const Eigen::Matrix<Scalar,ni,1> & ubi,
+    const Eigen::Matrix<Scalar,n,1> & lb,
+    const Eigen::Matrix<Scalar,n,1> & ub);
+  // min_x ½ xᵀ H x + xᵀf
+  // subject to:
+  //   A x = b
+  //   lb ≤ x ≤ u
+  //
+  // Templates:
+  //   Scalar  (e.g., double)
+  //   n  #H or Eigen::Dynamic if not known at compile time
+  //   m  #A or Eigen::Dynamic if not known at compile time
+  // Inputs:
+  //   H  #H by #H quadratic coefficients (only lower triangle used)
+  //   f  #H linear coefficients
+  //   A  #A by #H list of linear equality constraint coefficients
+  //   b  #A list of linear equality lower bounds
+  //   lb  #H list of lower bounds
+  //   ub  #H list of lower bounds
+  // Returns #H-long solution x
+  template <typename Scalar, int n, int m>
+  IGL_INLINE Eigen::Matrix<Scalar,n,1> quadprog(
+    const Eigen::Matrix<Scalar,n,n> & H,
+    const Eigen::Matrix<Scalar,n,1> & f,
+    const Eigen::Matrix<Scalar,m,n> & A,
+    const Eigen::Matrix<Scalar,m,1> & b,
+    const Eigen::Matrix<Scalar,n,1> & lb,
+    const Eigen::Matrix<Scalar,n,1> & ub);
+  // min_x ½ xᵀ H x + xᵀf
+  // subject to:
+  //   lb ≤ x ≤ u
+  //
+  // Templates:
+  //   Scalar  (e.g., double)
+  //   n  #H or Eigen::Dynamic if not known at compile time
+  // Inputs:
+  //   H  #H by #H quadratic coefficients (only lower triangle used)
+  //   f  #H linear coefficients
+  //   lb  #H list of lower bounds
+  //   ub  #H list of lower bounds
+  // Returns #H-long solution x
+  template <typename Scalar, int n>
+  IGL_INLINE Eigen::Matrix<Scalar,n,1> quadprog(
+    const Eigen::Matrix<Scalar,n,n> & H,
+    const Eigen::Matrix<Scalar,n,1> & f,
+    const Eigen::Matrix<Scalar,n,1> & lb,
+    const Eigen::Matrix<Scalar,n,1> & ub);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "quadprog.cpp"
+#endif 
+
+#endif 

include/igl/rigid_alignment.cpp

@@ -7,7 +7,6 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "rigid_alignment.h"
 #include "polar_svd.h"
-#include "matlab_format.h"
 #include <Eigen/Sparse>
 #include <Eigen/Cholesky>
 #include <vector>

include/igl/slice.cpp

@@ -189,6 +189,14 @@ IGL_INLINE DerivedX igl::slice(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::slice<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::DenseBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+// generated by autoexplicit.sh
+template void igl::slice<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::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, Eigen::Matrix<int, -1, -1, 0, -1, -1>&);
+// generated by autoexplicit.sh
+template void igl::slice<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Block<Eigen::Matrix<int, -1, -1, 0, -1, -1>, -1, 1, true>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::DenseBase<Eigen::Block<Eigen::Matrix<int, -1, -1, 0, -1, -1>, -1, 1, true> > const&, int, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
+// generated by autoexplicit.sh
+template void igl::slice<Eigen::VectorBlock<Eigen::Matrix<double, -1, 1, 0, -1, 1> const, -1>, Eigen::Block<Eigen::Matrix<int, -1, -1, 0, -1, -1>, -1, 1, true>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::VectorBlock<Eigen::Matrix<double, -1, 1, 0, -1, 1> const, -1> const&, Eigen::DenseBase<Eigen::Block<Eigen::Matrix<int, -1, -1, 0, -1, -1>, -1, 1, true> > const&, int, Eigen::Matrix<double, -1, 1, 0, -1, 1>&);
+// generated by autoexplicit.sh
 template void igl::slice<Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::DenseBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> > const&, Eigen::DenseBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template Eigen::Matrix<double, -1, -1, 0, -1, -1> igl::slice<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>>(Eigen::DenseBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const &, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>> const &, int);
 template Eigen::Matrix<int, -1, -1, 0, -1, -1> igl::slice<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>>(Eigen::DenseBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>> const &, Eigen::DenseBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>> const &);

include/igl/slice_mask.cpp

@@ -66,7 +66,11 @@ IGL_INLINE void igl::slice_mask(
         {
           if(R(i))
           {
-            Y.row(yi++) = X.row(i);
+            for(int j = 0;j<X.cols();j++)
+            {
+              Y(yi,j) = X(i,j);
+            }
+            yi++;
           }
         }
       }
@@ -84,7 +88,11 @@ IGL_INLINE void igl::slice_mask(
         {
           if(C(j))
           {
-            Y.col(yj++) = X.col(j);
+            for(int i = 0;i<X.rows();i++)
+            {
+              Y(i,yj) = X(i,j);
+            }
+            yj++;
           }
         }
       }
@@ -149,15 +157,16 @@ IGL_INLINE void igl::slice_mask(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
-// generated by autoexplicit.sh
+template void igl::slice_mask<Eigen::Matrix<double, 2, 1, 0, 2, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::DenseBase<Eigen::Matrix<double, 2, 1, 0, 2, 1> > const&, Eigen::Array<bool, -1, 1, 0, -1, 1> const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template void igl::slice_mask<Eigen::Matrix<double, 2, 2, 0, 2, 2>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::DenseBase<Eigen::Matrix<double, 2, 2, 0, 2, 2> > const&, Eigen::Array<bool, -1, 1, 0, -1, 1> const&, Eigen::Array<bool, -1, 1, 0, -1, 1> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::slice_mask<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::DenseBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::Array<bool, -1, 1, 0, -1, 1> const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+template void igl::slice_mask<Eigen::Matrix<double, 3, 3, 0, 3, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::DenseBase<Eigen::Matrix<double, 3, 3, 0, 3, 3> > const&, Eigen::Array<bool, -1, 1, 0, -1, 1> const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template void igl::slice_mask<Eigen::Matrix<double, 3, 1, 0, 3, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::DenseBase<Eigen::Matrix<double, 3, 1, 0, 3, 1> > const&, Eigen::Array<bool, -1, 1, 0, -1, 1> const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template void igl::slice_mask<Eigen::Matrix<double, 3, 3, 0, 3, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::DenseBase<Eigen::Matrix<double, 3, 3, 0, 3, 3> > const&, Eigen::Array<bool, -1, 1, 0, -1, 1> const&, Eigen::Array<bool, -1, 1, 0, -1, 1> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template Eigen::Array<int, -1, -1, 0, -1, -1> igl::slice_mask<Eigen::Array<int, -1, -1, 0, -1, -1> >(Eigen::DenseBase<Eigen::Array<int, -1, -1, 0, -1, -1> > const&, Eigen::Array<bool, -1, 1, 0, -1, 1> const&, int);
-// generated by autoexplicit.sh
 template Eigen::Matrix<int, -1, 3, 1, -1, 3> igl::slice_mask<Eigen::Matrix<int, -1, 3, 1, -1, 3> >(Eigen::DenseBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, Eigen::Array<bool, -1, 1, 0, -1, 1> const&, int);
-// generated by autoexplicit.sh
 template Eigen::Matrix<int, -1, 1, 0, -1, 1> igl::slice_mask<Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::Array<bool, -1, 1, 0, -1, 1> const&, int);
-// generated by autoexplicit.sh
 template Eigen::Array<int, -1, 3, 1, -1, 3> igl::slice_mask<Eigen::Array<int, -1, 3, 1, -1, 3> >(Eigen::DenseBase<Eigen::Array<int, -1, 3, 1, -1, 3> > const&, Eigen::Array<bool, -1, 1, 0, -1, 1> const&, int);
-// generated by autoexplicit.sh
 template void igl::slice_mask<bool, bool>(Eigen::SparseMatrix<bool, 0, int> const&, Eigen::Array<bool, -1, 1, 0, -1, 1> const&, int, Eigen::SparseMatrix<bool, 0, int>&);
 template void igl::slice_mask<Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<double, -1, 2, 0, -1, 2> >(Eigen::DenseBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> > const&, Eigen::Array<bool, -1, 1, 0, -1, 1> const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&);
 template void igl::slice_mask<Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 2, 0, -1, 2> >(Eigen::DenseBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, Eigen::Array<bool, -1, 1, 0, -1, 1> const&, int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&);

include/igl/sort.cpp

@@ -328,6 +328,8 @@ if(!ascending)
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::sort<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::DenseBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+// generated by autoexplicit.sh
 template void igl::sort<Eigen::Matrix<unsigned int, -1, 2, 0, -1, 2>, Eigen::Matrix<unsigned int, -1, -1, 0, -1, -1>, Eigen::Matrix<unsigned int, -1, -1, 0, -1, -1> >(Eigen::DenseBase<Eigen::Matrix<unsigned int, -1, 2, 0, -1, 2> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 0, -1, -1> >&);
 // generated by autoexplicit.sh
 template void igl::sort<Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::DenseBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);

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

include/igl/unique_edge_map.cpp

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

include/igl/voxel_grid.cpp

@@ -89,6 +89,12 @@ IGL_INLINE void igl::voxel_grid(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::voxel_grid<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 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, -1, -1, 0, -1, -1>::Scalar, int, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, 1, 3, 1, 1, 3> >&);
+// generated by autoexplicit.sh
+template void igl::voxel_grid<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<int, 1, 3, 1, 1, 3> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::Matrix<float, -1, 3, 1, -1, 3>::Scalar, int, int, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, 1, 3, 1, 1, 3> >&);
+// generated by autoexplicit.sh
+template void igl::voxel_grid<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, 1, 3, 1, 1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::Matrix<double, -1, 3, 1, -1, 3>::Scalar, int, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, 1, 3, 1, 1, 3> >&);
+// generated by autoexplicit.sh
 template void igl::voxel_grid<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, 1, 3, 1, 1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, int, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, 1, 3, 1, 1, 3> >&);
 // generated by autoexplicit.sh
 template void igl::voxel_grid<float, Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<int, 3, 1, 0, 3, 1> >(Eigen::AlignedBox<float, 3> const&, int, int, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, 3, 1, 0, 3, 1> >&);

include/igl/writePLY.cpp

@@ -413,6 +413,7 @@ template bool igl::writePLY<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matri
 template bool igl::writePLY<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, igl::FileEncoding);
 // generated by autoexplicit.sh
 template bool igl::writePLY<Eigen::Matrix<double, 8, 3, 0, 8, 3>, Eigen::Matrix<int, 12, 3, 0, 12, 3> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 8, 3, 0, 8, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, 12, 3, 0, 12, 3> > const&, igl::FileEncoding);
+template bool igl::writePLY<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&);
 template bool igl::writePLY<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&);
 template bool igl::writePLY<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, 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<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > 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<double, -1, -1, 0, -1, -1> > const&, std::vector<std::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::allocator<std::basic_string<char, std::char_traits<char>, std::allocator<char> > > > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, std::vector<std::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::allocator<std::basic_string<char, std::char_traits<char>, std::allocator<char> > > > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, std::vector<std::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::allocator<std::basic_string<char, std::char_traits<char>, std::allocator<char> > > > const&, std::vector<std::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::allocator<std::basic_string<char, std::char_traits<char>, std::allocator<char> > > > const&, igl::FileEncoding);
 template bool igl::writePLY<Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, std::vector<std::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::allocator<std::basic_string<char, std::char_traits<char>, std::allocator<char> > > > const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, std::vector<std::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::allocator<std::basic_string<char, std::char_traits<char>, std::allocator<char> > > > const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, std::vector<std::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::allocator<std::basic_string<char, std::char_traits<char>, std::allocator<char> > > > const&, std::vector<std::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::allocator<std::basic_string<char, std::char_traits<char>, std::allocator<char> > > > const&, igl::FileEncoding);

tests/include/igl/min_quad_with_fixed.cpp

@@ -0,0 +1,18 @@
+#include <test_common.h>
+#include <igl/min_quad_with_fixed.h>
+#include <igl/EPS.h>
+
+TEST_CASE("min_quad_with_fixed: dense", "[igl]" )
+{
+  const Eigen::Matrix<double,3,3> H = (Eigen::Matrix<double,3,3>(3,3)<<62,43,76,43,69,62,76,62,101).finished();
+  const Eigen::Matrix<double,3,1> f = (Eigen::Matrix<double,3,1>(3,1)<<9,8,5).finished();
+  const Eigen::Matrix<double,1,3> A = (Eigen::Matrix<double,1,3>(1,3)<<1,1,1).finished();
+  const Eigen::Matrix<double,1,1> b = (Eigen::Matrix<double,1,1>(1,1)<<2).finished();
+  const Eigen::Array<bool,3,1> k = (Eigen::Array<bool,3,1>()<<true,false,false).finished();
+  const Eigen::Matrix<double,3,1> bc = (Eigen::Matrix<double,3,1>(3,1)<<1,0,0).finished();
+  // Windows needs template args spelled out
+  const Eigen::Matrix<double,3,1> x = igl::min_quad_with_fixed<double,3,1>(H,f,k,bc,A,b);
+  REQUIRE(abs(x(0)- 1.0)<igl::EPS<double>());
+  REQUIRE(abs(x(1)- 1.5)<igl::EPS<double>());
+  REQUIRE(abs(x(2)- -.5)<igl::EPS<double>());
+}

tests/include/igl/polygon_corners.cpp

@@ -0,0 +1,28 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2020 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 <test_common.h>
+#include <igl/polygon_corners.h>
+#include <igl/matrix_to_list.h>
+#include <igl/matlab_format.h>
+#include <iostream>
+
+TEST_CASE("polygon_corners: quads", "[igl]")
+{
+  const Eigen::MatrixXi Q = (Eigen::MatrixXi(2,4)<< 0,1,2,3, 1,4,5,2).finished();
+  std::vector<std::vector<int>> vQ;
+  igl::matrix_to_list(Q,vQ);
+  const Eigen::VectorXi Iexact = (Eigen::VectorXi(8)<<0,1,2,3,1,4,5,2).finished();
+  const Eigen::VectorXi Cexact = (Eigen::VectorXi(3)<<0,4,8).finished();
+  Eigen::VectorXi I,C;
+  igl::polygon_corners(vQ,I,C);
+  test_common::assert_eq(I,Iexact);
+  test_common::assert_eq(C,Cexact);
+  igl::polygon_corners( Q,I,C);
+  test_common::assert_eq(I,Iexact);
+  test_common::assert_eq(C,Cexact);
+}

tests/include/igl/quadprog.cpp

@@ -0,0 +1,137 @@
+#include <test_common.h>
+#include <igl/quadprog.h>
+#include <igl/EPS.h>
+
+
+TEST_CASE("quadprog: box3", "[igl]" )
+{
+  {
+    Eigen::Matrix3d H = (Eigen::Matrix3d(3,3)<<
+        0.240548455386281,   0.237314308102107,  0.0436993831501944,
+        0.237314308102107,   0.326254049041854,-0.00021896091952631,
+        0.0436993831501944,-0.00021896091952631,   0.171175681280756
+        ).finished();
+    Eigen::Vector3d f = (Eigen::Vector3d(3,1)<<
+        0.222182612270718,
+        0.503254616893693,
+        -0.184619987497072
+        ).finished();
+    Eigen::Vector3d lb = (Eigen::Vector3d(3,1)<<
+        -1,
+        -1,
+        -1
+        ).finished();
+    Eigen::Vector3d ub = (Eigen::Vector3d(3,1)<<
+        1,
+        1,
+        1
+        ).finished(); 
+	// Windows needs template args spelled out
+    Eigen::Vector3d x = igl::quadprog<double,3>(H,f,lb,ub);
+    //std::cout<<igl::matlab_format(x,"x")<<std::endl;
+    REQUIRE(abs(x(0)- -0.118760635036839)<1e-7);
+    REQUIRE(abs(x(1)- -1)<1e-7);
+    REQUIRE(abs(x(2)- +1)<1e-7);
+  }
+  {
+    Eigen::Matrix3d H;
+    H<<1,0,0,0,1,0,0,0,1;
+    Eigen::Vector3d f( 0.5,-0.5,-0.5);
+    Eigen::Vector3d lb(0,0,0);
+    Eigen::Vector3d ub(1,1,1);
+	// Windows needs template args spelled out
+    Eigen::Vector3d x = igl::quadprog<double,3>(H,f,lb,ub);
+    REQUIRE(x(0)==0.0);
+    REQUIRE(x(1)==0.5);
+    REQUIRE(x(1)==0.5);
+  }
+  {
+    Eigen::Matrix3d H = (Eigen::Matrix3d(3,3)<<
+        1.06020279605748, 0.387347953430924,-0.653587847224834,
+        0.387347953430924, 0.323001970642631, -0.80259721932688,
+        -0.653587847224834, -0.80259721932688,  2.73709523329989
+        ).finished();
+    Eigen::Vector3d f = (Eigen::Vector3d(3,1)<<
+        -0.0155804986732503,
+        -0.00161174173383921,
+        -0.00903647945917485
+        ).finished();
+    Eigen::Vector3d lb = (Eigen::Vector3d(3,1)<<
+        0,
+        0,
+        0
+        ).finished();
+    Eigen::Vector3d ub = (Eigen::Vector3d(3,1)<<
+        0.015625,
+        0.015625,
+        0.015625
+        ).finished();
+	// Windows needs template args spelled out
+    Eigen::Vector3d x = igl::quadprog<double,3>(H,f,lb,ub);
+    Eigen::Vector3d xexact(0.015625, 0.013732474124087, 0.0110593284260843);
+    REQUIRE((x-xexact).array().abs().maxCoeff() < 1e-4);
+  }
+}
+
+TEST_CASE("quadprog: box2", "[igl]" )
+{
+  {
+    Eigen::Matrix2d H = (Eigen::Matrix2d(2,2)<<
+        0.683698654982294,-0.0521997092763332,
+        -0.0521997092763332,  0.800535063458999
+        ).finished();
+    Eigen::Vector2d f = (Eigen::Vector2d(2,1)<<
+        -0.733716332234936,
+        2.56401312736278
+        ).finished();
+    Eigen::Vector2d lb = (Eigen::Vector2d(2,1)<<
+        -1,
+        -1
+        ).finished();
+    Eigen::Vector2d ub = (Eigen::Vector2d(2,1)<<
+        1,
+        1
+        ).finished();
+	// Windows needs template args spelled out
+    Eigen::Vector2d x = igl::quadprog<double,2>(H,f,lb,ub);
+    //std::cout<<igl::matlab_format(x,"x")<<std::endl;
+    REQUIRE(abs(x(0)-0.99680848864073357)<1e-7);
+    REQUIRE(abs(x(1)- -1.)<1e-7);
+  }
+  {
+    Eigen::Matrix2d H = (Eigen::Matrix2d(2,2)<<
+        0.0708025678527926,-0.158030756288795,
+        -0.158030756288795, 0.360468620664163
+        ).finished();
+    Eigen::Vector2d f = (Eigen::Vector2d(2,1)<<
+        0.207229875768196,
+        -0.547595351878845
+        ).finished();
+    Eigen::Vector2d lb = (Eigen::Vector2d(2,1)<<
+        -1,
+        -1
+        ).finished();
+    Eigen::Vector2d ub = (Eigen::Vector2d(2,1)<<
+        1,
+        1
+        ).finished();
+	// Windows needs template args spelled out
+    Eigen::Vector2d x = igl::quadprog<double,2>(H,f,lb,ub);
+    //std::cout<<igl::matlab_format(x,"x")<<std::endl;
+    REQUIRE(abs(x(0)- -0.69487761491492939)<1e-7);
+    REQUIRE(abs(x(1)-1.0)<1e-7);
+  }
+  {
+    Eigen::Matrix2d H;
+    H<<0.4000,-0.2000,-0.2000,1.0000;
+    Eigen::Vector2d f(-0.3000,4.0000);
+    Eigen::Vector2d lb(0,0);
+    Eigen::Vector2d ub(1,1);
+	// Windows needs template args spelled out
+    Eigen::Vector2d x = igl::quadprog<double,2>(H,f,lb,ub);
+    //std::cout<<igl::matlab_format(x,"x")<<std::endl;
+    REQUIRE(abs(x(0)-0.75)<2e-16);
+    REQUIRE(abs(x(1)-0.0)<2e-16);
+  }
+}
+

tests/include/igl/sparse_voxel_grid.cpp

@@ -0,0 +1,22 @@
+#include <test_common.h>
+#include <igl/sparse_voxel_grid.h>
+#include <igl/unique_rows.h>
+
+TEST_CASE("sparse_voxel_grid: unique", "[igl]" )
+{
+  const std::function<double(const Eigen::RowVector3d & x)> f = 
+    [&](const Eigen::RowVector3d & x)->double
+  {
+    return x.norm() - 1.0;
+  };
+  Eigen::RowVector3d p0(0,1.0,0);
+  Eigen::MatrixXd GV;
+  Eigen::VectorXd Gf;
+  Eigen::Matrix<int,Eigen::Dynamic,8> GI;
+  igl::sparse_voxel_grid(p0,f,1,1024,Gf,GV,GI);
+  Eigen::MatrixXd uGV;
+  Eigen::VectorXi _1,_2;
+  igl::unique_rows(GV,uGV,_1,_2);
+  REQUIRE(GV.rows() == uGV.rows());
+}
+

tests/test_common.h

@@ -24,6 +24,28 @@
 #define IGL_DEBUG_OFF "[!hide]"
 #endif
 
+
+#include <igl/STR.h>
+template<>
+struct Catch::StringMaker<std::tuple<int,int,int> > 
+{
+  static std::string convert(std::tuple<int,int,int> const& t)
+  {
+    return 
+      STR("("<<std::get<0>(t)<<","<<std::get<1>(t)<<","<<std::get<2>(t)<<")");
+  }
+};
+template<>
+struct Catch::StringMaker<std::tuple<int,int,double> > 
+{
+  static std::string convert(std::tuple<int,int,double> const& t)
+  {
+    return 
+      STR("("<<std::get<0>(t)<<","<<std::get<1>(t)<<","<<std::get<2>(t)<<")");
+  }
+};
+
+
 namespace test_common
 {
   template<typename Param, typename Fun>

tutorial/715_MeshImplicitFunction/main.cpp

@@ -1,52 +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& x) -> double {
-    const double R = x.norm();
-    const double s = atan2(x.head(2).norm(),x(2));
-    const double p = atan2(x(1),x(0));
-    return pow(sin(s),2.)*(pow(cos(12.*s),3.)*0.1+pow(sin(6.*p),2)*0.2)+(R-0.5);
-    //return R-0.5;
+  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");
+
+  // 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();
   };
 
-  // We know that the point (0, 0, 0.5) lies on the implicit surface
-  Eigen::RowVector3d p0(0., 0., 0.5);
-
-  // Construct a sparse voxel grid whose cubes have edge length eps
-  // The cubes will form a thin shell around the implicit surface
-  const double eps = 0.01;
-
-  // CS will hold one scalar value at each cube vertex corresponding
-  // the value of the implicit at that vertex
-  Eigen::VectorXd CS;
-  // CV will hold the positions of the corners of the sparse voxel grid
-  Eigen::MatrixXd CV;
-  // 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;
-  // Construct the voxel grid, populating CS, CV, and CI
-  igl::sparse_voxel_grid(p0, scalar_func, eps, 200000/*expected_number_of_cubes*/, CS, CV, CI);
-
-  // Given the sparse voxel grid, use Marching Cubes to construct a triangle mesh of the surface
-  Eigen::MatrixXi F;
   Eigen::MatrixXd V;
-  igl::marching_cubes(CS, CV, CI, 0, 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();
+  Eigen::MatrixXi Q,F;
+  Eigen::MatrixXd mcV,mcN;
+  Eigen::MatrixXi mcF;
+
+  // 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)));
+
+  // 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
+    {
+
+      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);
+    }
+  }
+
+  // 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);
+  }
+
+  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();
 }
+