Merge branch 'master' of github.com:libigl/libigl into alecjacobson

Former-commit-id: 7be3e297d3d95706b0026824af30d3d3e38a4caf

include/igl/ConjugateFFSolverData.h

@@ -10,7 +10,9 @@
 #include "igl_inline.h"
 #include <Eigen/Core>
 #include <Eigen/Sparse>
-
+#include <igl/matlab_format.h>
+#include <iostream>
+using namespace std;
 namespace igl 
 {
   // Data class for the Conjugate Frame Field Solver
@@ -32,7 +34,7 @@ namespace igl
       Eigen::VectorXi indInteriorToFull;
       Eigen::VectorXi indFullToInterior;
 
-      Eigen::PlainObjectBase<DerivedV> B1, B2, FN;
+      DerivedV B1, B2, FN;
 
 
       Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic,1> kmin, kmax;
@@ -122,6 +124,8 @@ IGL_INLINE void igl::ConjugateFFSolverData<DerivedV, DerivedF>::computeCurvature
   kmax = kmax.bottomRows(numF);
   kmin = kmin.bottomRows(numF);
 
+  cerr<<igl::matlab_format(kmax,"kmax")<<endl;
+  cerr<<igl::matlab_format(kmin,"kmin")<<endl;
   //  kmax = dmax3.rowwise().norm();
   //  kmin = dmin3.rowwise().norm();
 

include/igl/colon.cpp

@@ -54,6 +54,7 @@ template void igl::colon<int, int, int>(int, int, Eigen::Matrix<int, -1, 1, 0, -
 template void igl::colon<int,long long int,int>(int,long long int,Eigen::Matrix<int,-1,1,0,-1,1> &);
 template void igl::colon<int, int, int, int>(int, int, int, Eigen::Matrix<int, -1, 1, 0, -1, 1>&);
 template void igl::colon<int, long, long>(int, long, Eigen::Matrix<long, -1, 1, 0, -1, 1>&);
+template void igl::colon<int, double, double, double>(int, double, double, Eigen::Matrix<double, -1, 1, 0, -1, 1>&);
 template void igl::colon<double, double, double>(double, double, Eigen::Matrix<double, -1, 1, 0, -1, 1>&);
 #ifdef WIN32
 template void igl::colon<int, long long,long>(int, long long, class Eigen::Matrix<long,-1,1,0,-1,1> &);

include/igl/comb_cross_field.cpp

@@ -27,13 +27,12 @@ namespace igl {
     const Eigen::PlainObjectBase<DerivedF> &F;
     const Eigen::PlainObjectBase<DerivedV> &PD1;
     const Eigen::PlainObjectBase<DerivedV> &PD2;
-//#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
-    Eigen::PlainObjectBase<DerivedV> N;
+    DerivedV N;
 
   private:
     // internal
-    Eigen::PlainObjectBase<DerivedF> TT;
-    Eigen::PlainObjectBase<DerivedF> TTi;
+    DerivedF TT;
+    DerivedF TTi;
 
 
   private:

include/igl/comb_line_field.cpp

@@ -25,14 +25,12 @@ public:
     const Eigen::PlainObjectBase<DerivedV> &V;
     const Eigen::PlainObjectBase<DerivedF> &F;
     const Eigen::PlainObjectBase<DerivedV> &PD1;
-//#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
-    Eigen::PlainObjectBase<DerivedV> N;
+    DerivedV N;
 
 private:
     // internal
-//#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
-    Eigen::PlainObjectBase<DerivedF> TT;
-    Eigen::PlainObjectBase<DerivedF> TTi;
+    DerivedF TT;
+    DerivedF TTi;
 
 
 private:

include/igl/conjugate_frame_fields.cpp

@@ -20,46 +20,45 @@ namespace igl {
   public:
     IGL_INLINE ConjugateFFSolver(const ConjugateFFSolverData<DerivedV, DerivedF> &_data,
                                  int _maxIter = 50,
-                                 const typename DerivedV::Scalar &_lambdaOrtho = .1,
-                                 const typename DerivedV::Scalar &_lambdaInit = 100,
-                                 const typename DerivedV::Scalar &_lambdaMultFactor = 1.01,
+                                 const typename DerivedV::Scalar _lambdaOrtho = .1,
+                                 const typename DerivedV::Scalar _lambdaInit = 100,
+                                 const typename DerivedV::Scalar _lambdaMultFactor = 1.01,
                                  bool _doHardConstraints = true);
-    IGL_INLINE bool solve(const Eigen::VectorXi &isConstrained,
-                          const Eigen::PlainObjectBase<DerivedO> &initialSolution,
-                          Eigen::PlainObjectBase<DerivedO> &output,
-                          typename DerivedV::Scalar *lambdaOut = NULL);
-
+    IGL_INLINE typename DerivedV::Scalar solve(const Eigen::VectorXi &isConstrained,
+                                               const Eigen::PlainObjectBase<DerivedO> &initialSolution,
+                                               Eigen::PlainObjectBase<DerivedO> &output);
+    
   private:
-
+    
     const ConjugateFFSolverData<DerivedV, DerivedF> &data;
-
+    
     //polyVF data
     Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> Acoeff, Bcoeff;
     Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 2> pvU, pvV;
     typename DerivedV::Scalar lambda;
-
+    
     //parameters
     typename DerivedV::Scalar lambdaOrtho;
     typename DerivedV::Scalar lambdaInit,lambdaMultFactor;
     int maxIter;
     bool doHardConstraints;
-
+    
     IGL_INLINE void localStep();
     IGL_INLINE void getPolyCoeffsForLocalSolve(const Eigen::Matrix<typename DerivedV::Scalar, 4, 1> &s,
                                                const Eigen::Matrix<typename DerivedV::Scalar, 4, 1> &z,
                                                Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 1> &polyCoeff);
-
+    
     IGL_INLINE void globalStep(const Eigen::Matrix<int, Eigen::Dynamic, 1>  &isConstrained,
                                const Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1>  &Ak,
                                const Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1>  &Bk);
     IGL_INLINE void minQuadWithKnownMini(const Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &Q,
-                         const Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &f,
-                         const Eigen::VectorXi isConstrained,
-                         const Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> &xknown,
+                                         const Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &f,
+                                         const Eigen::VectorXi isConstrained,
+                                         const Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> &xknown,
                                          Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> &x);
     IGL_INLINE void setFieldFromCoefficients();
     IGL_INLINE void setCoefficientsFromField();
-
+    
   };
 }
 
@@ -69,9 +68,9 @@ template <typename DerivedV, typename DerivedF, typename DerivedO>
 IGL_INLINE igl::ConjugateFFSolver<DerivedV, DerivedF, DerivedO>::
 ConjugateFFSolver(const ConjugateFFSolverData<DerivedV, DerivedF> &_data,
                   int _maxIter,
-                  const typename DerivedV::Scalar &_lambdaOrtho,
-                  const typename DerivedV::Scalar &_lambdaInit,
-                  const typename DerivedV::Scalar &_lambdaMultFactor,
+                  const typename DerivedV::Scalar _lambdaOrtho,
+                  const typename DerivedV::Scalar _lambdaInit,
+                  const typename DerivedV::Scalar _lambdaMultFactor,
                   bool _doHardConstraints):
 data(_data),
 lambdaOrtho(_lambdaOrtho),
@@ -102,7 +101,7 @@ getPolyCoeffsForLocalSolve(const Eigen::Matrix<typename DerivedV::Scalar, 4, 1>
   typename DerivedV::Scalar z1 = z(1);
   typename DerivedV::Scalar z2 = z(2);
   typename DerivedV::Scalar z3 = z(3);
-
+  
   polyCoeff.resize(7,1);
   polyCoeff(0) =  s0*s0* s1*s1* s2*s2* s3* z3*z3 +  s0*s0* s1*s1* s2* s3*s3* z2*z2 +  s0*s0* s1* s2*s2* s3*s3* z1*z1 +  s0* s1*s1* s2*s2* s3*s3* z0*z0 ;
   polyCoeff(1) = 2* s0*s0* s1*s1* s2* s3* z2*z2 + 2* s0*s0* s1*s1* s2* s3* z3*z3 + 2* s0*s0* s1* s2*s2* s3* z1*z1 + 2* s0*s0* s1* s2*s2* s3* z3*z3 + 2* s0*s0* s1* s2* s3*s3* z1*z1 + 2* s0*s0* s1* s2* s3*s3* z2*z2 + 2* s0* s1*s1* s2*s2* s3* z0*z0 + 2* s0* s1*s1* s2*s2* s3* z3*z3 + 2* s0* s1*s1* s2* s3*s3* z0*z0 + 2* s0* s1*s1* s2* s3*s3* z2*z2 + 2* s0* s1* s2*s2* s3*s3* z0*z0 + 2* s0* s1* s2*s2* s3*s3* z1*z1 ;
@@ -111,7 +110,7 @@ getPolyCoeffsForLocalSolve(const Eigen::Matrix<typename DerivedV::Scalar, 4, 1>
   polyCoeff(4) =  s0*s0* s1* z1*z1 +  s0*s0* s2* z2*z2 +  s0*s0* s3* z3*z3 +  s0* s1*s1* z0*z0 + 4* s0* s1* s2* z0*z0 + 4* s0* s1* s2* z1*z1 + 4* s0* s1* s2* z2*z2 + 4* s0* s1* s3* z0*z0 + 4* s0* s1* s3* z1*z1 + 4* s0* s1* s3* z3*z3 +  s0* s2*s2* z0*z0 + 4* s0* s2* s3* z0*z0 + 4* s0* s2* s3* z2*z2 + 4* s0* s2* s3* z3*z3 +  s0* s3*s3* z0*z0 +  s1*s1* s2* z2*z2 +  s1*s1* s3* z3*z3 +  s1* s2*s2* z1*z1 + 4* s1* s2* s3* z1*z1 + 4* s1* s2* s3* z2*z2 + 4* s1* s2* s3* z3*z3 +  s1* s3*s3* z1*z1 +  s2*s2* s3* z3*z3 +  s2* s3*s3* z2*z2;
   polyCoeff(5) = 2* s0* s1* z0*z0 + 2* s0* s1* z1*z1 + 2* s0* s2* z0*z0 + 2* s0* s2* z2*z2 + 2* s0* s3* z0*z0 + 2* s0* s3* z3*z3 + 2* s1* s2* z1*z1 + 2* s1* s2* z2*z2 + 2* s1* s3* z1*z1 + 2* s1* s3* z3*z3 + 2* s2* s3* z2*z2 + 2* s2* s3* z3*z3 ;
   polyCoeff(6) =  s0* z0*z0 +  s1* z1*z1 +  s2* z2*z2 +  s3* z3*z3;
-
+  
 }
 
 
@@ -124,12 +123,12 @@ localStep()
     Eigen::Matrix<typename DerivedV::Scalar, 4, 1> xproj; xproj << pvU.row(j).transpose(),pvV.row(j).transpose();
     Eigen::Matrix<typename DerivedV::Scalar, 4, 1> z = data.UH[j].transpose()*xproj;
     Eigen::Matrix<typename DerivedV::Scalar, 4, 1> x;
-
+    
     Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 1> polyCoeff;
     getPolyCoeffsForLocalSolve(data.s[j], z, polyCoeff);
     Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> roots;
     igl::polyRoots<typename DerivedV::Scalar, typename DerivedV::Scalar> (polyCoeff, roots );
-
+    
     //  find closest real root to xproj
     typename DerivedV::Scalar minDist = 1e10;
     for (int i =0; i< 6; ++i)
@@ -144,9 +143,9 @@ localStep()
         minDist = dist;
         x = candidate;
       }
-
+      
     }
-
+    
     pvU.row(j) << x(0),x(1);
     pvV.row(j) << x(2),x(3);
   }
@@ -174,15 +173,15 @@ globalStep(const Eigen::Matrix<int, Eigen::Dynamic, 1>  &isConstrained,
            const Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1>  &Bk)
 {
   setCoefficientsFromField();
-
+  
   Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > I;
   igl::speye(data.numF, data.numF, I);
   Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > QA = data.DDA+lambda*data.planarityWeight+lambdaOrtho*I;
   Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > fA = (-2*lambda*data.planarityWeight*Acoeff).sparseView();
-
+  
   Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > QB = data.DDB+lambda*data.planarityWeight;
   Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > fB = (-2*lambda*data.planarityWeight*Bcoeff).sparseView();
-
+  
   if(doHardConstraints)
   {
     minQuadWithKnownMini(QA, fA, isConstrained, Ak, Acoeff);
@@ -196,7 +195,7 @@ globalStep(const Eigen::Matrix<int, Eigen::Dynamic, 1>  &isConstrained,
     minQuadWithKnownMini(QB, fB, isknown_, xknown_, Bcoeff);
   }
   setFieldFromCoefficients();
-
+  
 }
 
 
@@ -210,10 +209,10 @@ setFieldFromCoefficients()
     //    matlab code from roots (given there are no trailing zeros in the polynomial coefficients)
     Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> polyCoeff(5,1);
     polyCoeff<<1., 0., -Acoeff(i), 0., Bcoeff(i);
-
+    
     Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> roots;
     polyRoots<std::complex<typename DerivedV::Scalar>>(polyCoeff,roots);
-
+    
     std::complex<typename DerivedV::Scalar> u = roots[0];
     int maxi = -1;
     float maxd = -1;
@@ -230,7 +229,7 @@ setFieldFromCoefficients()
     pvU(i,0) = real(u); pvU(i,1) = imag(u);
     pvV(i,0) = real(v); pvV(i,1) = imag(v);
   }
-
+  
 }
 
 template<typename DerivedV, typename DerivedF, typename DerivedO>
@@ -242,11 +241,11 @@ minQuadWithKnownMini(const Eigen::SparseMatrix<std::complex<typename DerivedV::S
                      Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> &x)
 {
   int N = Q.rows();
-
+  
   int nc = xknown.rows();
   Eigen::VectorXi known; known.setZero(nc,1);
   Eigen::VectorXi unknown; unknown.setZero(N-nc,1);
-
+  
   int indk = 0, indu = 0;
   for (int i = 0; i<N; ++i)
     if (isConstrained[i])
@@ -259,21 +258,21 @@ minQuadWithKnownMini(const Eigen::SparseMatrix<std::complex<typename DerivedV::S
       unknown[indu] = i;
       indu++;
     }
-
+  
   Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar>> Quu, Quk;
-
+  
   igl::slice(Q,unknown, unknown, Quu);
   igl::slice(Q,unknown, known, Quk);
-
-
+  
+  
   std::vector<typename Eigen::Triplet<std::complex<typename DerivedV::Scalar> > > tripletList;
-
+  
   Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > fu(N-nc,1);
-
+  
   igl::slice(f,unknown, Eigen::VectorXi::Zero(1,1), fu);
-
+  
   Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > rhs = (Quk*xknown).sparseView()+.5*fu;
-
+  
   Eigen::SparseLU< Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar>>> solver;
   solver.compute(-Quu);
   if(solver.info()!=Eigen::Success)
@@ -287,7 +286,7 @@ minQuadWithKnownMini(const Eigen::SparseMatrix<std::complex<typename DerivedV::S
     std::cerr<<"Solving failed!"<<std::endl;
     return;
   }
-
+  
   indk = 0, indu = 0;
   x.setZero(N,1);
   for (int i = 0; i<N; ++i)
@@ -295,21 +294,20 @@ minQuadWithKnownMini(const Eigen::SparseMatrix<std::complex<typename DerivedV::S
       x[i] = xknown[indk++];
     else
       x[i] = b.coeff(indu++,0);
-
+  
 }
 
 
 template<typename DerivedV, typename DerivedF, typename DerivedO>
-IGL_INLINE bool igl::ConjugateFFSolver<DerivedV, DerivedF, DerivedO>::
+IGL_INLINE typename DerivedV::Scalar igl::ConjugateFFSolver<DerivedV, DerivedF, DerivedO>::
 solve(const Eigen::VectorXi &isConstrained,
       const Eigen::PlainObjectBase<DerivedO> &initialSolution,
-      Eigen::PlainObjectBase<DerivedO> &output,
-      typename DerivedV::Scalar *lambdaOut)
+      Eigen::PlainObjectBase<DerivedO> &output)
 {
   int numConstrained = isConstrained.sum();
   // coefficient values
   Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> Ak, Bk;
-
+  
   pvU.resize(data.numF,2);
   pvV.resize(data.numF,2);
   for (int fi = 0; fi <data.numF; ++fi)
@@ -334,54 +332,54 @@ solve(const Eigen::VectorXi &isConstrained,
       ind ++;
     }
   }
-
-
-
+  
+  
+  
   typename DerivedV::Scalar smoothnessValue;
   Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 1> conjValues;
   typename DerivedV::Scalar meanConj;
   typename DerivedV::Scalar maxConj;
-
+  
   data.evaluateConjugacy(pvU, pvV, conjValues);
   meanConj = conjValues.cwiseAbs().mean();
   maxConj = conjValues.cwiseAbs().maxCoeff();
   printf("Initial max non-conjugacy: %.5g\n",maxConj);
-
+  
   smoothnessValue = (Acoeff.adjoint()*data.DDA*Acoeff + Bcoeff.adjoint()*data.DDB*Bcoeff).real()[0];
   printf("\n\nInitial smoothness: %.5g\n",smoothnessValue);
-
+  
   lambda = lambdaInit;
-
+  
   bool doit = false;
   for (int iter = 0; iter<maxIter; ++iter)
   {
     printf("\n\n--- Iteration %d ---\n",iter);
-
+    
     typename DerivedV::Scalar oldMeanConj = meanConj;
-
+    
     localStep();
     globalStep(isConstrained, Ak, Bk);
-
-
+    
+    
     smoothnessValue = (Acoeff.adjoint()*data.DDA*Acoeff + Bcoeff.adjoint()*data.DDB*Bcoeff).real()[0];
-
+    
     printf("Smoothness: %.5g\n",smoothnessValue);
-
+    
     data.evaluateConjugacy(pvU, pvV, conjValues);
     meanConj = conjValues.cwiseAbs().mean();
     maxConj = conjValues.cwiseAbs().maxCoeff();
     printf("Mean/Max non-conjugacy: %.5g, %.5g\n",meanConj,maxConj);
     typename DerivedV::Scalar diffMeanConj = fabs(oldMeanConj-meanConj);
-
+    
     if (diffMeanConj<1e-4)
       doit = true;
-
+    
     if (doit)
       lambda = lambda*lambdaMultFactor;
     printf(" %d %.5g %.5g\n",iter, smoothnessValue,maxConj);
-
+    
   }
-
+  
   output.setZero(data.numF,6);
   for (int fi=0; fi<data.numF; ++fi)
   {
@@ -390,26 +388,23 @@ solve(const Eigen::VectorXi &isConstrained,
     output.block(fi,0, 1, 3) = pvU(fi,0)*b1 + pvU(fi,1)*b2;
     output.block(fi,3, 1, 3) = pvV(fi,0)*b1 + pvV(fi,1)*b2;
   }
-
-  if (lambdaOut)
-    *lambdaOut = lambda;
-
-  return true;
+  
+  return lambda;
 }
 
 
 
 template <typename DerivedV, typename DerivedF, typename DerivedO>
 IGL_INLINE void igl::conjugate_frame_fields(const Eigen::PlainObjectBase<DerivedV> &V,
-                                            const Eigen::PlainObjectBase<DerivedF> &F,
-                                            const Eigen::VectorXi &isConstrained,
-                                            const Eigen::PlainObjectBase<DerivedO> &initialSolution,
-                                            Eigen::PlainObjectBase<DerivedO> &output,
-                                            int maxIter,
-                                            const typename DerivedV::Scalar &lambdaOrtho,
-                                            const typename DerivedV::Scalar &lambdaInit,
-                                            const typename DerivedV::Scalar &lambdaMultFactor,
-                                            bool doHardConstraints)
+                                       const Eigen::PlainObjectBase<DerivedF> &F,
+                                       const Eigen::VectorXi &isConstrained,
+                                       const Eigen::PlainObjectBase<DerivedO> &initialSolution,
+                                       Eigen::PlainObjectBase<DerivedO> &output,
+                                       int maxIter,
+                                       const typename DerivedV::Scalar lambdaOrtho,
+                                       const typename DerivedV::Scalar lambdaInit,
+                                       const typename DerivedV::Scalar lambdaMultFactor,
+                                       bool doHardConstraints)
 {
   igl::ConjugateFFSolverData<DerivedV, DerivedF> csdata(V, F);
   igl::ConjugateFFSolver<DerivedV, DerivedF, DerivedO> cs(csdata, maxIter, lambdaOrtho, lambdaInit, lambdaMultFactor, doHardConstraints);
@@ -417,22 +412,22 @@ IGL_INLINE void igl::conjugate_frame_fields(const Eigen::PlainObjectBase<Derived
 }
 
 template <typename DerivedV, typename DerivedF, typename DerivedO>
-IGL_INLINE void igl::conjugate_frame_fields(const igl::ConjugateFFSolverData<DerivedV, DerivedF> &csdata,
-                                            const Eigen::VectorXi &isConstrained,
-                                            const Eigen::PlainObjectBase<DerivedO> &initialSolution,
-                                            Eigen::PlainObjectBase<DerivedO> &output,
-                                            int maxIter,
-                                            const typename DerivedV::Scalar &lambdaOrtho,
-                                            const typename DerivedV::Scalar &lambdaInit,
-                                            const typename DerivedV::Scalar &lambdaMultFactor,
-                                            bool doHardConstraints,
-                                            typename DerivedV::Scalar *lambdaOut)
+IGL_INLINE typename DerivedV::Scalar igl::conjugate_frame_fields(const igl::ConjugateFFSolverData<DerivedV, DerivedF> &csdata,
+                                                                 const Eigen::VectorXi &isConstrained,
+                                                                 const Eigen::PlainObjectBase<DerivedO> &initialSolution,
+                                                                 Eigen::PlainObjectBase<DerivedO> &output,
+                                                                 int maxIter,
+                                                                 const typename DerivedV::Scalar lambdaOrtho,
+                                                                 const typename DerivedV::Scalar lambdaInit,
+                                                                 const typename DerivedV::Scalar lambdaMultFactor,
+                                                                 bool doHardConstraints)
 {
   igl::ConjugateFFSolver<DerivedV, DerivedF, DerivedO> cs(csdata, maxIter, lambdaOrtho, lambdaInit, lambdaMultFactor, doHardConstraints);
-  cs.solve(isConstrained, initialSolution, output, lambdaOut);
+  typename DerivedV::Scalar lambdaOut = cs.solve(isConstrained, initialSolution, output);
+  return lambdaOut;
 }
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-template void igl::conjugate_frame_fields<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(igl::ConjugateFFSolverData<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, int, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar const&, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar const&, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar const&, bool, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar*);
+template Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar igl::conjugate_frame_fields<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(igl::ConjugateFFSolverData<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, int, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, bool);
 #endif

include/igl/conjugate_frame_fields.h

@@ -17,17 +17,7 @@ namespace igl {
   //todo
   // TODO: isConstrained should become a list of indices for consistency with
   //       n_polyvector
-  /// Given 2 vectors centered on origin calculate the rotation matrix from first to the second
-
-  // Inputs:
-  //   v0, v1         the two #3 by 1 vectors
-  //   normalized     boolean, if false, then the vectors are normalized prior to the calculation
-  // Output:
-  //                  3 by 3 rotation matrix that takes v0 to v1
-  //
-  template <typename DerivedV, typename DerivedF>
-  class ConjugateFFSolverData;
-
+  
   template <typename DerivedV, typename DerivedF, typename DerivedO>
   IGL_INLINE void conjugate_frame_fields(const Eigen::PlainObjectBase<DerivedV> &V,
                                          const Eigen::PlainObjectBase<DerivedF> &F,
@@ -35,23 +25,22 @@ namespace igl {
                                          const Eigen::PlainObjectBase<DerivedO> &initialSolution,
                                          Eigen::PlainObjectBase<DerivedO> &output,
                                          int _maxIter = 50,
-                                         const typename DerivedV::Scalar &_lambdaOrtho = .1,
-                                         const typename DerivedV::Scalar &_lambdaInit = 100,
-                                         const typename DerivedV::Scalar &_lambdaMultFactor = 1.01,
+                                         const typename DerivedV::Scalar _lambdaOrtho = .1,
+                                         const typename DerivedV::Scalar _lambdaInit = 100,
+                                         const typename DerivedV::Scalar _lambdaMultFactor = 1.01,
                                          bool _doHardConstraints = true);
-
+  
   template <typename DerivedV, typename DerivedF, typename DerivedO>
-  IGL_INLINE void conjugate_frame_fields(const ConjugateFFSolverData<DerivedV, DerivedF> &csdata,
-                                         const Eigen::VectorXi &isConstrained,
-                                         const Eigen::PlainObjectBase<DerivedO> &initialSolution,
-                                         Eigen::PlainObjectBase<DerivedO> &output,
-                                         int _maxIter = 50,
-                                         const typename DerivedV::Scalar &_lambdaOrtho = .1,
-                                         const typename DerivedV::Scalar &_lambdaInit = 100,
-                                         const typename DerivedV::Scalar &_lambdaMultFactor = 1.01,
-                                         bool _doHardConstraints = true,
-                                         typename DerivedV::Scalar *lambdaOut = NULL);
-
+  IGL_INLINE typename DerivedV::Scalar conjugate_frame_fields(const ConjugateFFSolverData<DerivedV, DerivedF> &csdata,
+                                                              const Eigen::VectorXi &isConstrained,
+                                                              const Eigen::PlainObjectBase<DerivedO> &initialSolution,
+                                                              Eigen::PlainObjectBase<DerivedO> &output,
+                                                              int _maxIter = 50,
+                                                              const typename DerivedV::Scalar _lambdaOrtho = .1,
+                                                              const typename DerivedV::Scalar _lambdaInit = 100,
+                                                              const typename DerivedV::Scalar _lambdaMultFactor = 1.01,
+                                                              bool _doHardConstraints = true);
+  
 };
 
 

include/igl/copyleft/comiso/miq.cpp

@@ -297,18 +297,18 @@ namespace comiso {
   private:
     const Eigen::PlainObjectBase<DerivedV> &V;
     const Eigen::PlainObjectBase<DerivedF> &F;
-    Eigen::PlainObjectBase<DerivedV> Vcut;
-    Eigen::PlainObjectBase<DerivedF> Fcut;
+    DerivedV Vcut;
+    DerivedF Fcut;
     Eigen::MatrixXd UV_out;
-    Eigen::PlainObjectBase<DerivedF> FUV_out;
+    DerivedF FUV_out;
 
     // internal
-    Eigen::PlainObjectBase<DerivedF> TT;
-    Eigen::PlainObjectBase<DerivedF> TTi;
+    DerivedF TT;
+    DerivedF TTi;
 
     // Stiffness per face
     Eigen::VectorXd Handle_Stiffness;
-    Eigen::PlainObjectBase<DerivedV> B1, B2, B3;
+    DerivedV B1, B2, B3;
 
   public:
     IGL_INLINE MIQ_class(const Eigen::PlainObjectBase<DerivedV> &V_,
@@ -1486,13 +1486,13 @@ IGL_INLINE void igl::copyleft::comiso::miq(
     std::vector<std::vector<int> > hardFeatures)
 {
 
-  Eigen::PlainObjectBase<DerivedV> BIS1, BIS2;
+  DerivedV BIS1, BIS2;
   igl::compute_frame_field_bisectors(V, F, PD1, PD2, BIS1, BIS2);
 
-  Eigen::PlainObjectBase<DerivedV> BIS1_combed, BIS2_combed;
+  DerivedV BIS1_combed, BIS2_combed;
   igl::comb_cross_field(V, F, BIS1, BIS2, BIS1_combed, BIS2_combed);
 
-  Eigen::PlainObjectBase<DerivedF> Handle_MMatch;
+  DerivedF Handle_MMatch;
   igl::cross_field_missmatch(V, F, BIS1_combed, BIS2_combed, true, Handle_MMatch);
 
   Eigen::Matrix<int, Eigen::Dynamic, 1> isSingularity, singularityIndex;
@@ -1501,7 +1501,7 @@ IGL_INLINE void igl::copyleft::comiso::miq(
   Eigen::Matrix<int, Eigen::Dynamic, 3> Handle_Seams;
   igl::cut_mesh_from_singularities(V, F, Handle_MMatch, Handle_Seams);
 
-  Eigen::PlainObjectBase<DerivedV> PD1_combed, PD2_combed;
+  DerivedV PD1_combed, PD2_combed;
   igl::comb_frame_field(V, F, PD1, PD2, BIS1_combed, BIS2_combed, PD1_combed, PD2_combed);
 
   igl::copyleft::comiso::miq(V,

include/igl/cross_field_missmatch.cpp

@@ -28,17 +28,17 @@ namespace igl {
     const Eigen::PlainObjectBase<DerivedF> &F;
     const Eigen::PlainObjectBase<DerivedV> &PD1;
     const Eigen::PlainObjectBase<DerivedV> &PD2;
-//#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
-    Eigen::PlainObjectBase<DerivedV> N;
+    
+    DerivedV N;
 
   private:
     // internal
     std::vector<bool> V_border; // bool
     std::vector<std::vector<int> > VF;
     std::vector<std::vector<int> > VFi;
-//#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
-    Eigen::PlainObjectBase<DerivedF> TT;
-    Eigen::PlainObjectBase<DerivedF> TTi;
+    
+    DerivedF TT;
+    DerivedF TTi;
 
 
   private:

include/igl/cut_mesh.cpp

@@ -29,8 +29,7 @@ namespace igl {
     int num_scalar_variables;
 
     // per face indexes of vertex in the solver
-//#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
-    Eigen::PlainObjectBase<DerivedF> HandleS_Index;
+    DerivedF HandleS_Index;
 
     // per vertex variable indexes
     std::vector<std::vector<int> > HandleV_Integer;

include/igl/cut_mesh_from_singularities.cpp

@@ -30,9 +30,8 @@ namespace igl {
     const Eigen::PlainObjectBase<DerivedM> &Handle_MMatch;
 
     Eigen::VectorXi F_visited;
-//#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
-    Eigen::PlainObjectBase<DerivedF> TT;
-    Eigen::PlainObjectBase<DerivedF> TTi;
+    DerivedF TT;
+    DerivedF TTi;
 
     Eigen::MatrixXi E, F2E, E2F;
   protected:

include/igl/facet_components.cpp

@@ -88,4 +88,5 @@ IGL_INLINE void igl::facet_components(
 // Explicit template specialization
 template void igl::facet_components<long, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<long, -1, 1, 0, -1, 1> >(std::vector<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >, std::allocator<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
 template void igl::facet_components<int, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(std::vector<std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >, std::allocator<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> >&);
+template void igl::facet_components<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/find.cpp

@@ -126,4 +126,5 @@ template void igl::find<double, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matr
 template void igl::find<double, 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> >&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 // generated by autoexplicit.sh
 template void igl::find<double, 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::SparseMatrix<double, 0, int> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::find<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/hsv_to_rgb.cpp

@@ -68,4 +68,5 @@ template void igl::hsv_to_rgb<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::M
 template void igl::hsv_to_rgb<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::hsv_to_rgb<Eigen::Matrix<unsigned char, 64, 3, 1, 64, 3>, Eigen::Matrix<unsigned char, 64, 3, 1, 64, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<unsigned char, 64, 3, 1, 64, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned char, 64, 3, 1, 64, 3> >&);
 template void igl::hsv_to_rgb<Eigen::Matrix<float, 64, 3, 1, 64, 3>, Eigen::Matrix<float, 64, 3, 1, 64, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<float, 64, 3, 1, 64, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, 64, 3, 1, 64, 3> >&);
+template void igl::hsv_to_rgb<double>(double const*, double*);
 #endif

include/igl/integrable_polyvector_fields.cpp

@@ -190,8 +190,7 @@ makeFieldCCW(Eigen::MatrixXd &sol3D)
   {
     //take all 4 vectors (including opposites) and pick two that are in ccw order
     all << sol3D.row(fi), -sol3D.row(fi);
-	Eigen::VectorXi inv_order_unused;
-    igl::sort_vectors_ccw(all, FN.row(fi).eval(), order, true, t, false, inv_order_unused);
+	  igl::sort_vectors_ccw(all, FN.row(fi).eval(), order, t);
     //if we are in a constrained face, we need to make sure that the first vector is always the same vector as in the constraints
     if(is_constrained_face[fi])
     {

include/igl/integrable_polyvector_fields.h

@@ -133,7 +133,7 @@ public:
   //per-edge angles (for parallel transport)
   Eigen::VectorXd K;
   //local bases
-  Eigen::PlainObjectBase<DerivedV> B1, B2, FN;
+  DerivedV B1, B2, FN;
 
   //Solver Data
   Eigen::VectorXd residuals;

include/igl/is_vertex_manifold.cpp

@@ -97,4 +97,5 @@ IGL_INLINE bool igl::is_vertex_manifold(
 
 #ifdef IGL_STATIC_LIBRARY
 template bool igl::is_vertex_manifold<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template bool igl::is_vertex_manifold<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/matlab_format.cpp

@@ -54,7 +54,7 @@ igl::matlab_format(
   if(!name.empty())
   {
     prefix = name + "IJV = ";
-    suffix = "\n"+name + " = sparse("+name+"IJV(:,1),"+name+"IJV(:,2),"+name+"IJV(:,3));";
+    suffix = "\n"+name + " = sparse("+name+"IJV(:,1),"+name+"IJV(:,2),"+name+"IJV(:,3),"+std::to_string(S.rows())+","+std::to_string(S.cols())+" );";
   }
   return STR(""<<
     SIJV.format(Eigen::IOFormat(
@@ -128,4 +128,5 @@ template Eigen::WithFormat<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const igl::m
 template Eigen::WithFormat<Eigen::Matrix<double, 2, 3, 0, 2, 3> > const igl::matlab_format<Eigen::Matrix<double, 2, 3, 0, 2, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 2, 3, 0, 2, 3> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
 template Eigen::WithFormat<Eigen::Matrix<double, 3, 2, 0, 3, 2> > const igl::matlab_format<Eigen::Matrix<double, 3, 2, 0, 3, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 2, 0, 3, 2> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
 template Eigen::WithFormat<Eigen::Matrix<float, -1, 1, 0, -1, 1> > const igl::matlab_format<Eigen::Matrix<float, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 1, 0, -1, 1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+template Eigen::WithFormat<Eigen::Matrix<int, 2, 2, 0, 2, 2> > const igl::matlab_format<Eigen::Matrix<int, 2, 2, 0, 2, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<int, 2, 2, 0, 2, 2> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
 #endif

include/igl/n_polyvector.cpp

@@ -39,8 +39,7 @@ namespace igl {
     Eigen::VectorXi indInteriorToFull;
     Eigen::VectorXi indFullToInterior;
 
-//#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
-    Eigen::PlainObjectBase<DerivedV> B1, B2, FN;
+    DerivedV B1, B2, FN;
 
     IGL_INLINE void computek();
     IGL_INLINE void setFieldFromGeneralCoefficients(const  std::vector<Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1> > &coeffs,

include/igl/n_polyvector_general.cpp

@@ -7,465 +7,462 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 
 #include <igl/n_polyvector_general.h>
-//#include <igl/edge_topology.h>
-//#include <igl/local_basis.h>
-//#include <igl/nchoosek.h>
-//#include <igl/slice.h>
-//#include <igl/polyroots.h>
-//#include <Eigen/Sparse>
-//#include <Eigen/Geometry>
-//#include <iostream>
-//#include <iostream>
-//
-//namespace igl {
-//  template <typename DerivedV, typename DerivedF>
-//  class GeneralPolyVectorFieldFinder
-//  {
-//  private:
-//    const Eigen::PlainObjectBase<DerivedV> &V;
-//    const Eigen::PlainObjectBase<DerivedF> &F; int numF;
-//    const int n;
-//
-//    Eigen::MatrixXi EV; int numE;
-//    Eigen::MatrixXi F2E;
-//    Eigen::MatrixXi E2F;
-//    Eigen::VectorXd K;
-//
-//    Eigen::VectorXi isBorderEdge;
-//    int numInteriorEdges;
-//    Eigen::Matrix<int,Eigen::Dynamic,2> E2F_int;
-//    Eigen::VectorXi indInteriorToFull;
-//    Eigen::VectorXi indFullToInterior;
-//
-////#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
-//    Eigen::PlainObjectBase<DerivedV> B1, B2, FN;
-//
-//    IGL_INLINE void computek();
-//    IGL_INLINE void setFieldFromGeneralCoefficients(const  std::vector<Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1>> &coeffs,
-//                                                    std::vector<Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 2> > &pv);
-//    IGL_INLINE void computeCoefficientLaplacian(int n, Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &D);
-//    IGL_INLINE void getGeneralCoeffConstraints(const Eigen::VectorXi &isConstrained,
-//                                    const Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &cfW,
-//                                    int k,
-//                                    const Eigen::VectorXi &rootsIndex,
-//                                    Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1> &Ck);
-//    IGL_INLINE void precomputeInteriorEdges();
-//
-//
-//    IGL_INLINE void minQuadWithKnownMini(const Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &Q,
-//                                         const Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &f,
-//                                         const Eigen::VectorXi isConstrained,
-//                                         const Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> &xknown,
-//                                         Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> &x);
-//
-//  public:
-//    IGL_INLINE GeneralPolyVectorFieldFinder(const Eigen::PlainObjectBase<DerivedV> &_V,
-//                                     const Eigen::PlainObjectBase<DerivedF> &_F,
-//                                     const int &_n);
-//    IGL_INLINE bool solve(const Eigen::VectorXi &isConstrained,
-//               const Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &cfW,
-//               const Eigen::VectorXi &rootsIndex,
-//               Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &output);
-//
-//  };
-//}
-//
-//template<typename DerivedV, typename DerivedF>
-//IGL_INLINE igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::
-//          GeneralPolyVectorFieldFinder(const Eigen::PlainObjectBase<DerivedV> &_V,
-//                                const Eigen::PlainObjectBase<DerivedF> &_F,
-//                                const int &_n):
-//V(_V),
-//F(_F),
-//numF(_F.rows()),
-//n(_n)
-//{
-//
-//  igl::edge_topology(V,F,EV,F2E,E2F);
-//  numE = EV.rows();
-//
-//
-//  precomputeInteriorEdges();
-//
-//  igl::local_basis(V,F,B1,B2,FN);
-//
-//  computek();
-//
-//};
-//
-//
-//template<typename DerivedV, typename DerivedF>
-//IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::
-//precomputeInteriorEdges()
-//{
-//  // Flag border edges
-//  numInteriorEdges = 0;
-//  isBorderEdge.setZero(numE,1);
-//  indFullToInterior = -1*Eigen::VectorXi::Ones(numE,1);
-//
-//  for(unsigned i=0; i<numE; ++i)
-//  {
-//    if ((E2F(i,0) == -1) || ((E2F(i,1) == -1)))
-//      isBorderEdge[i] = 1;
-//      else
-//      {
-//        indFullToInterior[i] = numInteriorEdges;
-//        numInteriorEdges++;
-//      }
-//  }
-//
-//  E2F_int.resize(numInteriorEdges, 2);
-//  indInteriorToFull.setZero(numInteriorEdges,1);
-//  int ii = 0;
-//  for (int k=0; k<numE; ++k)
-//  {
-//    if (isBorderEdge[k])
-//      continue;
-//    E2F_int.row(ii) = E2F.row(k);
-//    indInteriorToFull[ii] = k;
-//    ii++;
-//  }
-//
-//}
-//
-//
-//template<typename DerivedV, typename DerivedF>
-//IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::
-//minQuadWithKnownMini(const Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &Q,
-//                          const Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &f,
-//                     const Eigen::VectorXi isConstrained,
-//                          const Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> &xknown,
-//                          Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> &x)
-//{
-//  int N = Q.rows();
-//
-//  int nc = xknown.rows();
-//  Eigen::VectorXi known; known.setZero(nc,1);
-//  Eigen::VectorXi unknown; unknown.setZero(N-nc,1);
-//
-//  int indk = 0, indu = 0;
-//  for (int i = 0; i<N; ++i)
-//    if (isConstrained[i])
-//    {
-//      known[indk] = i;
-//      indk++;
-//    }
-//    else
-//    {
-//      unknown[indu] = i;
-//      indu++;
-//    }
-//
-//  Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar>> Quu, Quk;
-//
-//  igl::slice(Q,unknown, unknown, Quu);
-//  igl::slice(Q,unknown, known, Quk);
-//
-//
-//  std::vector<typename Eigen::Triplet<std::complex<typename DerivedV::Scalar> > > tripletList;
-//
-//  Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > fu(N-nc,1);
-//
-//  igl::slice(f,unknown, Eigen::VectorXi::Zero(1,1), fu);
-//
-//  Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > rhs = (Quk*xknown).sparseView()+.5*fu;
-//
-//  Eigen::SparseLU< Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar>>> solver;
-//  solver.compute(-Quu);
-//  if(solver.info()!=Eigen::Success)
-//  {
-//    std::cerr<<"Decomposition failed!"<<std::endl;
-//    return;
-//  }
-//  Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar>>  b  = solver.solve(rhs);
-//  if(solver.info()!=Eigen::Success)
-//  {
-//    std::cerr<<"Solving failed!"<<std::endl;
-//    return;
-//  }
-//
-//  indk = 0, indu = 0;
-//  x.setZero(N,1);
-//  for (int i = 0; i<N; ++i)
-//    if (isConstrained[i])
-//      x[i] = xknown[indk++];
-//    else
-//      x[i] = b.coeff(indu++,0);
-//
-//}
-//
-//
-//
-//template<typename DerivedV, typename DerivedF>
-//IGL_INLINE bool igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::
-//                     solve(const Eigen::VectorXi &isConstrained,
-//                           const Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &cfW,
-//                           const Eigen::VectorXi &rootsIndex,
-//                           Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &output)
-//{
-//	std::cerr << "This code is broken!" << std::endl;
-//	exit(1);
-//  // polynomial is of the form:
-//  // z^(2n) +
-//  // -c[0]z^(2n-1) +
-//  // c[1]z^(2n-2) +
-//  // -c[2]z^(2n-3) +
-//  // ... +
-//  // (-1)^n c[n-1]
-//  //std::vector<Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1>> coeffs(n,Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1>::Zero(numF, 1));
-//  //for (int i =0; i<n; ++i)
-//  //{
-//  //  int degree = i+1;
-//  //  Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1> Ck;
-//  //  getGeneralCoeffConstraints(isConstrained,
-//  //                             cfW,
-//  //                             i,
-//  //                             rootsIndex,
-//  //                             Ck);
-//
-//  //  Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > DD;
-//  //  computeCoefficientLaplacian(degree, DD);
-//  //  Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > f; f.resize(numF,1);
-//
-//  //  if (isConstrained.sum() == numF)
-//  //    coeffs[i] = Ck;
-//  //  else
-//  //    minQuadWithKnownMini(DD, f, isConstrained, Ck, coeffs[i]);
-//  //}
-//  //std::vector<Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 2> > pv;
-//  //setFieldFromGeneralCoefficients(coeffs, pv);
-//  //output.setZero(numF,3*n);
-//  //for (int fi=0; fi<numF; ++fi)
-//  //{
-//  //  const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> &b1 = B1.row(fi);
-//  //  const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> &b2 = B2.row(fi);
-//  //  for (int i=0; i<n; ++i)
-//  //    output.block(fi,3*i, 1, 3) = pv[i](fi,0)*b1 + pv[i](fi,1)*b2;
-//  //}
-//  return true;
-//}
-//
-//template<typename DerivedV, typename DerivedF>
-//IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::setFieldFromGeneralCoefficients(const  std::vector<Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1>> &coeffs,
-//                                                            std::vector<Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 2>> &pv)
-//{
-//  pv.assign(n, Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 2>::Zero(numF, 2));
-//  for (int i = 0; i <numF; ++i)
-//  {
-//
-//    //    poly coefficients: 1, 0, -Acoeff, 0, Bcoeff
-//    //    matlab code from roots (given there are no trailing zeros in the polynomial coefficients)
-//    Eigen::Matrix<std::complex<DerivedV::Scalar>, Eigen::Dynamic,1> polyCoeff;
-//    polyCoeff.setZero(n+1,1);
-//    polyCoeff[0] = 1.;
-//    int sign = 1;
-//    for (int k =0; k<n; ++k)
-//    {
-//      sign = -sign;
-//      int degree = k+1;
-//      polyCoeff[degree] = (1.*sign)*coeffs[k](i);
-//    }
-//
-//    Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1> roots;
-//    igl::polyRoots<std::complex<typename DerivedV::Scalar>, typename DerivedV::Scalar >(polyCoeff,roots);
-//    for (int k=0; k<n; ++k)
-//    {
-//      pv[k](i,0) = real(roots[k]);
-//      pv[k](i,1) = imag(roots[k]);
-//    }
-//  }
-//
-//}
-//
-//
-//template<typename DerivedV, typename DerivedF>
-//IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::computeCoefficientLaplacian(int n, Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &D)
-//{
-//  std::vector<Eigen::Triplet<std::complex<typename DerivedV::Scalar> >> tripletList;
-//
-//  // For every non-border edge
-//  for (unsigned eid=0; eid<numE; ++eid)
-//  {
-//    if (!isBorderEdge[eid])
-//    {
-//      int fid0 = E2F(eid,0);
-//      int fid1 = E2F(eid,1);
-//
-//      tripletList.push_back(Eigen::Triplet<std::complex<typename DerivedV::Scalar> >(fid0,
-//                                           fid0,
-//                                           std::complex<typename DerivedV::Scalar>(1.)));
-//      tripletList.push_back(Eigen::Triplet<std::complex<typename DerivedV::Scalar> >(fid1,
-//                                           fid1,
-//                                           std::complex<typename DerivedV::Scalar>(1.)));
-//      tripletList.push_back(Eigen::Triplet<std::complex<typename DerivedV::Scalar> >(fid0,
-//                                           fid1,
-//                                                                                     -1.*std::polar(1.,-1.*n*K[eid])));
-//      tripletList.push_back(Eigen::Triplet<std::complex<typename DerivedV::Scalar> >(fid1,
-//                                           fid0,
-//                                                                                     -1.*std::polar(1.,1.*n*K[eid])));
-//
-//    }
-//  }
-//  D.resize(numF,numF);
-//  D.setFromTriplets(tripletList.begin(), tripletList.end());
-//
-//
-//}
-//
-////this gives the coefficients without the (-1)^k that multiplies them
-//template<typename DerivedV, typename DerivedF>
-//IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::getGeneralCoeffConstraints(const Eigen::VectorXi &isConstrained,
-//                                                       const Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &cfW,
-//                                                       int k,
-//                                                       const Eigen::VectorXi &rootsIndex,
-//                                                       Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1> &Ck)
-//{
-//  int numConstrained = isConstrained.sum();
-//  Ck.resize(numConstrained,1);
-//  // int n = rootsIndex.cols();
-//
-//  Eigen::MatrixXi allCombs;
-//  {
-//    Eigen::VectorXi V = Eigen::VectorXi::LinSpaced(n,0,n-1);
-//    igl::nchoosek(V,k+1,allCombs);
-//  }
-//
-//  int ind = 0;
-//  for (int fi = 0; fi <numF; ++fi)
-//  {
-//    const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> &b1 = B1.row(fi);
-//    const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> &b2 = B2.row(fi);
-//    if(isConstrained[fi])
-//    {
-//      std::complex<typename DerivedV::Scalar> ck(0);
-//
-//      for (int j = 0; j < allCombs.rows(); ++j)
-//      {
-//        std::complex<typename DerivedV::Scalar> tk(1.);
-//        //collect products
-//        for (int i = 0; i < allCombs.cols(); ++i)
-//        {
-//          int index = allCombs(j,i);
-//
-//          int ri = rootsIndex[index];
-//          Eigen::Matrix<typename DerivedV::Scalar, 1, 3> w;
-//          if (ri>0)
-//            w = cfW.block(fi,3*(ri-1),1,3);
-//          else
-//            w = -cfW.block(fi,3*(-ri-1),1,3);
-//          typename DerivedV::Scalar w0 = w.dot(b1);
-//          typename DerivedV::Scalar w1 = w.dot(b2);
-//          std::complex<typename DerivedV::Scalar> u(w0,w1);
-//          tk*= u;
-//        }
-//        //collect sum
-//        ck += tk;
-//      }
-//      Ck(ind) = ck;
-//      ind ++;
-//    }
-//  }
-//
-//
-//}
-//
-//template<typename DerivedV, typename DerivedF>
-//IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::computek()
-//{
-//  K.setZero(numE);
-//  // For every non-border edge
-//  for (unsigned eid=0; eid<numE; ++eid)
-//  {
-//    if (!isBorderEdge[eid])
-//    {
-//      int fid0 = E2F(eid,0);
-//      int fid1 = E2F(eid,1);
-//
-//      Eigen::Matrix<typename DerivedV::Scalar, 1, 3> N0 = FN.row(fid0);
-//      Eigen::Matrix<typename DerivedV::Scalar, 1, 3> N1 = FN.row(fid1);
-//
-//      // find common edge on triangle 0 and 1
-//      int fid0_vc = -1;
-//      int fid1_vc = -1;
-//      for (unsigned i=0;i<3;++i)
-//      {
-//        if (F2E(fid0,i) == eid)
-//          fid0_vc = i;
-//        if (F2E(fid1,i) == eid)
-//          fid1_vc = i;
-//      }
-//      assert(fid0_vc != -1);
-//      assert(fid1_vc != -1);
-//
-//      Eigen::Matrix<typename DerivedV::Scalar, 1, 3> common_edge = V.row(F(fid0,(fid0_vc+1)%3)) - V.row(F(fid0,fid0_vc));
-//      common_edge.normalize();
-//
-//      // Map the two triangles in a new space where the common edge is the x axis and the N0 the z axis
-//      Eigen::Matrix<typename DerivedV::Scalar, 3, 3> P;
-//      Eigen::Matrix<typename DerivedV::Scalar, 1, 3> o = V.row(F(fid0,fid0_vc));
-//      Eigen::Matrix<typename DerivedV::Scalar, 1, 3> tmp = -N0.cross(common_edge);
-//      P << common_edge, tmp, N0;
-////      P.transposeInPlace();
-//
-//
-//      Eigen::Matrix<typename DerivedV::Scalar, 3, 3> V0;
-//      V0.row(0) = V.row(F(fid0,0)) -o;
-//      V0.row(1) = V.row(F(fid0,1)) -o;
-//      V0.row(2) = V.row(F(fid0,2)) -o;
-//
-//      V0 = (P*V0.transpose()).transpose();
-//
-////      assert(V0(0,2) < 1e-10);
-////      assert(V0(1,2) < 1e-10);
-////      assert(V0(2,2) < 1e-10);
-//
-//      Eigen::Matrix<typename DerivedV::Scalar, 3, 3> V1;
-//      V1.row(0) = V.row(F(fid1,0)) -o;
-//      V1.row(1) = V.row(F(fid1,1)) -o;
-//      V1.row(2) = V.row(F(fid1,2)) -o;
-//      V1 = (P*V1.transpose()).transpose();
-//
-////      assert(V1(fid1_vc,2) < 10e-10);
-////      assert(V1((fid1_vc+1)%3,2) < 10e-10);
-//
-//      // compute rotation R such that R * N1 = N0
-//      // i.e. map both triangles to the same plane
-//      double alpha = -atan2(V1((fid1_vc+2)%3,2),V1((fid1_vc+2)%3,1));
-//
-//      Eigen::Matrix<typename DerivedV::Scalar, 3, 3> R;
-//      R << 1,          0,            0,
-//      0, cos(alpha), -sin(alpha) ,
-//      0, sin(alpha),  cos(alpha);
-//      V1 = (R*V1.transpose()).transpose();
-//
-////      assert(V1(0,2) < 1e-10);
-////      assert(V1(1,2) < 1e-10);
-////      assert(V1(2,2) < 1e-10);
-//
-//      // measure the angle between the reference frames
-//      // k_ij is the angle between the triangle on the left and the one on the right
-//      Eigen::Matrix<typename DerivedV::Scalar, 1, 3> ref0 = V0.row(1) - V0.row(0);
-//      Eigen::Matrix<typename DerivedV::Scalar, 1, 3> ref1 = V1.row(1) - V1.row(0);
-//
-//      ref0.normalize();
-//      ref1.normalize();
-//
-//      double ktemp = atan2(ref1(1),ref1(0)) - atan2(ref0(1),ref0(0));
-//
-//      // just to be sure, rotate ref0 using angle ktemp...
-//      Eigen::Matrix<typename DerivedV::Scalar, 2, 2> R2;
-//      R2 << cos(ktemp), -sin(ktemp), sin(ktemp), cos(ktemp);
-//
-//      Eigen::Matrix<typename DerivedV::Scalar, 1, 2> tmp1 = R2*(ref0.head(2)).transpose();
-//
-////      assert(tmp1(0) - ref1(0) < 1e-10);
-////      assert(tmp1(1) - ref1(1) < 1e-10);
-//
-//      K[eid] = ktemp;
-//    }
-//  }
-//
-//}
+#include <igl/edge_topology.h>
+#include <igl/local_basis.h>
+#include <igl/nchoosek.h>
+#include <igl/slice.h>
+#include <igl/polyroots.h>
+#include <Eigen/Sparse>
+#include <Eigen/Geometry>
+#include <iostream>
+#include <iostream>
+
+namespace igl {
+  template <typename DerivedV, typename DerivedF>
+  class GeneralPolyVectorFieldFinder
+  {
+  private:
+    const Eigen::PlainObjectBase<DerivedV> &V;
+    const Eigen::PlainObjectBase<DerivedF> &F; int numF;
+    const int n;
+
+    Eigen::MatrixXi EV; int numE;
+    Eigen::MatrixXi F2E;
+    Eigen::MatrixXi E2F;
+    Eigen::VectorXd K;
+
+    Eigen::VectorXi isBorderEdge;
+    int numInteriorEdges;
+    Eigen::Matrix<int,Eigen::Dynamic,2> E2F_int;
+    Eigen::VectorXi indInteriorToFull;
+    Eigen::VectorXi indFullToInterior;
+
+    DerivedV B1, B2, FN;
+
+    IGL_INLINE void computek();
+    IGL_INLINE void setFieldFromGeneralCoefficients(const  std::vector<Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1>> &coeffs,
+                                                    std::vector<Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 2> > &pv);
+    IGL_INLINE void computeCoefficientLaplacian(int n, Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &D);
+    IGL_INLINE void getGeneralCoeffConstraints(const Eigen::VectorXi &isConstrained,
+                                    const Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &cfW,
+                                    int k,
+                                    const Eigen::VectorXi &rootsIndex,
+                                    Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1> &Ck);
+    IGL_INLINE void precomputeInteriorEdges();
+
+
+    IGL_INLINE void minQuadWithKnownMini(const Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &Q,
+                                         const Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &f,
+                                         const Eigen::VectorXi isConstrained,
+                                         const Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> &xknown,
+                                         Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> &x);
+
+  public:
+    IGL_INLINE GeneralPolyVectorFieldFinder(const Eigen::PlainObjectBase<DerivedV> &_V,
+                                     const Eigen::PlainObjectBase<DerivedF> &_F,
+                                     const int &_n);
+    IGL_INLINE bool solve(const Eigen::VectorXi &isConstrained,
+               const Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &cfW,
+               const Eigen::VectorXi &rootsIndex,
+               Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &output);
+
+  };
+}
+
+template<typename DerivedV, typename DerivedF>
+IGL_INLINE igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::
+          GeneralPolyVectorFieldFinder(const Eigen::PlainObjectBase<DerivedV> &_V,
+                                const Eigen::PlainObjectBase<DerivedF> &_F,
+                                const int &_n):
+V(_V),
+F(_F),
+numF(_F.rows()),
+n(_n)
+{
+
+  igl::edge_topology(V,F,EV,F2E,E2F);
+  numE = EV.rows();
+
+
+  precomputeInteriorEdges();
+
+  igl::local_basis(V,F,B1,B2,FN);
+
+  computek();
+
+};
+
+
+template<typename DerivedV, typename DerivedF>
+IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::
+precomputeInteriorEdges()
+{
+  // Flag border edges
+  numInteriorEdges = 0;
+  isBorderEdge.setZero(numE,1);
+  indFullToInterior = -1*Eigen::VectorXi::Ones(numE,1);
+
+  for(unsigned i=0; i<numE; ++i)
+  {
+    if ((E2F(i,0) == -1) || ((E2F(i,1) == -1)))
+      isBorderEdge[i] = 1;
+      else
+      {
+        indFullToInterior[i] = numInteriorEdges;
+        numInteriorEdges++;
+      }
+  }
+
+  E2F_int.resize(numInteriorEdges, 2);
+  indInteriorToFull.setZero(numInteriorEdges,1);
+  int ii = 0;
+  for (int k=0; k<numE; ++k)
+  {
+    if (isBorderEdge[k])
+      continue;
+    E2F_int.row(ii) = E2F.row(k);
+    indInteriorToFull[ii] = k;
+    ii++;
+  }
+
+}
+
+
+template<typename DerivedV, typename DerivedF>
+IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::
+minQuadWithKnownMini(const Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &Q,
+                          const Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &f,
+                     const Eigen::VectorXi isConstrained,
+                          const Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> &xknown,
+                          Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> &x)
+{
+  int N = Q.rows();
+
+  int nc = xknown.rows();
+  Eigen::VectorXi known; known.setZero(nc,1);
+  Eigen::VectorXi unknown; unknown.setZero(N-nc,1);
+
+  int indk = 0, indu = 0;
+  for (int i = 0; i<N; ++i)
+    if (isConstrained[i])
+    {
+      known[indk] = i;
+      indk++;
+    }
+    else
+    {
+      unknown[indu] = i;
+      indu++;
+    }
+
+  Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar>> Quu, Quk;
+
+  igl::slice(Q,unknown, unknown, Quu);
+  igl::slice(Q,unknown, known, Quk);
+
+
+  std::vector<typename Eigen::Triplet<std::complex<typename DerivedV::Scalar> > > tripletList;
+
+  Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > fu(N-nc,1);
+
+  igl::slice(f,unknown, Eigen::VectorXi::Zero(1,1), fu);
+
+  Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > rhs = (Quk*xknown).sparseView()+.5*fu;
+
+  Eigen::SparseLU< Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar>>> solver;
+  solver.compute(-Quu);
+  if(solver.info()!=Eigen::Success)
+  {
+    std::cerr<<"Decomposition failed!"<<std::endl;
+    return;
+  }
+  Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar>>  b  = solver.solve(rhs);
+  if(solver.info()!=Eigen::Success)
+  {
+    std::cerr<<"Solving failed!"<<std::endl;
+    return;
+  }
+
+  indk = 0, indu = 0;
+  x.setZero(N,1);
+  for (int i = 0; i<N; ++i)
+    if (isConstrained[i])
+      x[i] = xknown[indk++];
+    else
+      x[i] = b.coeff(indu++,0);
+
+}
+
+
+
+template<typename DerivedV, typename DerivedF>
+IGL_INLINE bool igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::
+                     solve(const Eigen::VectorXi &isConstrained,
+                           const Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &cfW,
+                           const Eigen::VectorXi &rootsIndex,
+                           Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &output)
+{
+  // polynomial is of the form:
+  // z^(2n) +
+  // -c[0]z^(2n-1) +
+  // c[1]z^(2n-2) +
+  // -c[2]z^(2n-3) +
+  // ... +
+  // (-1)^n c[n-1]
+  std::vector<Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1>> coeffs(n,Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1>::Zero(numF, 1));
+  for (int i =0; i<n; ++i)
+  {
+    int degree = i+1;
+    Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1> Ck;
+    getGeneralCoeffConstraints(isConstrained,
+                               cfW,
+                               i,
+                               rootsIndex,
+                               Ck);
+
+    Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > DD;
+    computeCoefficientLaplacian(degree, DD);
+    Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > f; f.resize(numF,1);
+
+    if (isConstrained.sum() == numF)
+      coeffs[i] = Ck;
+    else
+      minQuadWithKnownMini(DD, f, isConstrained, Ck, coeffs[i]);
+  }
+  std::vector<Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 2> > pv;
+  setFieldFromGeneralCoefficients(coeffs, pv);
+  output.setZero(numF,3*n);
+  for (int fi=0; fi<numF; ++fi)
+  {
+    const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> &b1 = B1.row(fi);
+    const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> &b2 = B2.row(fi);
+    for (int i=0; i<n; ++i)
+      output.block(fi,3*i, 1, 3) = pv[i](fi,0)*b1 + pv[i](fi,1)*b2;
+  }
+  return true;
+}
+
+template<typename DerivedV, typename DerivedF>
+IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::setFieldFromGeneralCoefficients(const  std::vector<Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1>> &coeffs,
+                                                            std::vector<Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 2>> &pv)
+{
+  pv.assign(n, Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 2>::Zero(numF, 2));
+  for (int i = 0; i <numF; ++i)
+  {
+
+    //    poly coefficients: 1, 0, -Acoeff, 0, Bcoeff
+    //    matlab code from roots (given there are no trailing zeros in the polynomial coefficients)
+    Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1> polyCoeff;
+    polyCoeff.setZero(n+1,1);
+    polyCoeff[0] = 1.;
+    int sign = 1;
+    for (int k =0; k<n; ++k)
+    {
+      sign = -sign;
+      int degree = k+1;
+      polyCoeff[degree] = (1.*sign)*coeffs[k](i);
+    }
+
+    Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1> roots;
+    igl::polyRoots<std::complex<typename DerivedV::Scalar>, typename DerivedV::Scalar >(polyCoeff,roots);
+    for (int k=0; k<n; ++k)
+    {
+      pv[k](i,0) = real(roots[k]);
+      pv[k](i,1) = imag(roots[k]);
+    }
+  }
+
+}
+
+
+template<typename DerivedV, typename DerivedF>
+IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::computeCoefficientLaplacian(int n, Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &D)
+{
+  std::vector<Eigen::Triplet<std::complex<typename DerivedV::Scalar> >> tripletList;
+
+  // For every non-border edge
+  for (unsigned eid=0; eid<numE; ++eid)
+  {
+    if (!isBorderEdge[eid])
+    {
+      int fid0 = E2F(eid,0);
+      int fid1 = E2F(eid,1);
+
+      tripletList.push_back(Eigen::Triplet<std::complex<typename DerivedV::Scalar> >(fid0,
+                                           fid0,
+                                           std::complex<typename DerivedV::Scalar>(1.)));
+      tripletList.push_back(Eigen::Triplet<std::complex<typename DerivedV::Scalar> >(fid1,
+                                           fid1,
+                                           std::complex<typename DerivedV::Scalar>(1.)));
+      tripletList.push_back(Eigen::Triplet<std::complex<typename DerivedV::Scalar> >(fid0,
+                                           fid1,
+                                                                                     -1.*std::polar(1.,-1.*n*K[eid])));
+      tripletList.push_back(Eigen::Triplet<std::complex<typename DerivedV::Scalar> >(fid1,
+                                           fid0,
+                                                                                     -1.*std::polar(1.,1.*n*K[eid])));
+
+    }
+  }
+  D.resize(numF,numF);
+  D.setFromTriplets(tripletList.begin(), tripletList.end());
+
+
+}
+
+//this gives the coefficients without the (-1)^k that multiplies them
+template<typename DerivedV, typename DerivedF>
+IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::getGeneralCoeffConstraints(const Eigen::VectorXi &isConstrained,
+                                                       const Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &cfW,
+                                                       int k,
+                                                       const Eigen::VectorXi &rootsIndex,
+                                                       Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1> &Ck)
+{
+  int numConstrained = isConstrained.sum();
+  Ck.resize(numConstrained,1);
+  // int n = rootsIndex.cols();
+
+  Eigen::MatrixXi allCombs;
+  {
+    Eigen::VectorXi V = Eigen::VectorXi::LinSpaced(n,0,n-1);
+    igl::nchoosek(V,k+1,allCombs);
+  }
+
+  int ind = 0;
+  for (int fi = 0; fi <numF; ++fi)
+  {
+    const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> &b1 = B1.row(fi);
+    const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> &b2 = B2.row(fi);
+    if(isConstrained[fi])
+    {
+      std::complex<typename DerivedV::Scalar> ck(0);
+
+      for (int j = 0; j < allCombs.rows(); ++j)
+      {
+        std::complex<typename DerivedV::Scalar> tk(1.);
+        //collect products
+        for (int i = 0; i < allCombs.cols(); ++i)
+        {
+          int index = allCombs(j,i);
+
+          int ri = rootsIndex[index];
+          Eigen::Matrix<typename DerivedV::Scalar, 1, 3> w;
+          if (ri>0)
+            w = cfW.block(fi,3*(ri-1),1,3);
+          else
+            w = -cfW.block(fi,3*(-ri-1),1,3);
+          typename DerivedV::Scalar w0 = w.dot(b1);
+          typename DerivedV::Scalar w1 = w.dot(b2);
+          std::complex<typename DerivedV::Scalar> u(w0,w1);
+          tk*= u;
+        }
+        //collect sum
+        ck += tk;
+      }
+      Ck(ind) = ck;
+      ind ++;
+    }
+  }
+
+
+}
+
+template<typename DerivedV, typename DerivedF>
+IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::computek()
+{
+  K.setZero(numE);
+  // For every non-border edge
+  for (unsigned eid=0; eid<numE; ++eid)
+  {
+    if (!isBorderEdge[eid])
+    {
+      int fid0 = E2F(eid,0);
+      int fid1 = E2F(eid,1);
+
+      Eigen::Matrix<typename DerivedV::Scalar, 1, 3> N0 = FN.row(fid0);
+      Eigen::Matrix<typename DerivedV::Scalar, 1, 3> N1 = FN.row(fid1);
+
+      // find common edge on triangle 0 and 1
+      int fid0_vc = -1;
+      int fid1_vc = -1;
+      for (unsigned i=0;i<3;++i)
+      {
+        if (F2E(fid0,i) == eid)
+          fid0_vc = i;
+        if (F2E(fid1,i) == eid)
+          fid1_vc = i;
+      }
+      assert(fid0_vc != -1);
+      assert(fid1_vc != -1);
+
+      Eigen::Matrix<typename DerivedV::Scalar, 1, 3> common_edge = V.row(F(fid0,(fid0_vc+1)%3)) - V.row(F(fid0,fid0_vc));
+      common_edge.normalize();
+
+      // Map the two triangles in a new space where the common edge is the x axis and the N0 the z axis
+      Eigen::Matrix<typename DerivedV::Scalar, 3, 3> P;
+      Eigen::Matrix<typename DerivedV::Scalar, 1, 3> o = V.row(F(fid0,fid0_vc));
+      Eigen::Matrix<typename DerivedV::Scalar, 1, 3> tmp = -N0.cross(common_edge);
+      P << common_edge, tmp, N0;
+//      P.transposeInPlace();
+
+
+      Eigen::Matrix<typename DerivedV::Scalar, 3, 3> V0;
+      V0.row(0) = V.row(F(fid0,0)) -o;
+      V0.row(1) = V.row(F(fid0,1)) -o;
+      V0.row(2) = V.row(F(fid0,2)) -o;
+
+      V0 = (P*V0.transpose()).transpose();
+
+//      assert(V0(0,2) < 1e-10);
+//      assert(V0(1,2) < 1e-10);
+//      assert(V0(2,2) < 1e-10);
+
+      Eigen::Matrix<typename DerivedV::Scalar, 3, 3> V1;
+      V1.row(0) = V.row(F(fid1,0)) -o;
+      V1.row(1) = V.row(F(fid1,1)) -o;
+      V1.row(2) = V.row(F(fid1,2)) -o;
+      V1 = (P*V1.transpose()).transpose();
+
+//      assert(V1(fid1_vc,2) < 10e-10);
+//      assert(V1((fid1_vc+1)%3,2) < 10e-10);
+
+      // compute rotation R such that R * N1 = N0
+      // i.e. map both triangles to the same plane
+      double alpha = -atan2(V1((fid1_vc+2)%3,2),V1((fid1_vc+2)%3,1));
+
+      Eigen::Matrix<typename DerivedV::Scalar, 3, 3> R;
+      R << 1,          0,            0,
+      0, cos(alpha), -sin(alpha) ,
+      0, sin(alpha),  cos(alpha);
+      V1 = (R*V1.transpose()).transpose();
+
+//      assert(V1(0,2) < 1e-10);
+//      assert(V1(1,2) < 1e-10);
+//      assert(V1(2,2) < 1e-10);
+
+      // measure the angle between the reference frames
+      // k_ij is the angle between the triangle on the left and the one on the right
+      Eigen::Matrix<typename DerivedV::Scalar, 1, 3> ref0 = V0.row(1) - V0.row(0);
+      Eigen::Matrix<typename DerivedV::Scalar, 1, 3> ref1 = V1.row(1) - V1.row(0);
+
+      ref0.normalize();
+      ref1.normalize();
+
+      double ktemp = atan2(ref1(1),ref1(0)) - atan2(ref0(1),ref0(0));
+
+      // just to be sure, rotate ref0 using angle ktemp...
+      Eigen::Matrix<typename DerivedV::Scalar, 2, 2> R2;
+      R2 << cos(ktemp), -sin(ktemp), sin(ktemp), cos(ktemp);
+
+      Eigen::Matrix<typename DerivedV::Scalar, 1, 2> tmp1 = R2*(ref0.head(2)).transpose();
+
+//      assert(tmp1(0) - ref1(0) < 1e-10);
+//      assert(tmp1(1) - ref1(1) < 1e-10);
+
+      K[eid] = ktemp;
+    }
+  }
+
+}
 
 
 IGL_INLINE void igl::n_polyvector_general(const Eigen::MatrixXd &V,
@@ -475,19 +472,17 @@ IGL_INLINE void igl::n_polyvector_general(const Eigen::MatrixXd &V,
                              const Eigen::VectorXi &I,
                              Eigen::MatrixXd &output)
 {
-	// This functions is broken, please contact Olga Diamanti
-	assert(0);
-  //Eigen::VectorXi isConstrained = Eigen::VectorXi::Constant(F.rows(),0);
-  //Eigen::MatrixXd cfW = Eigen::MatrixXd::Constant(F.rows(),bc.cols(),0);
-
-  //for(unsigned i=0; i<b.size();++i)
-  //{
-  //  isConstrained(b(i)) = 1;
-  //  cfW.row(b(i)) << bc.row(i);
-  //}
-  //int n = I.rows();
-  //igl::GeneralPolyVectorFieldFinder<Eigen::MatrixXd, Eigen::MatrixXi> pvff(V,F,n);
-  //pvff.solve(isConstrained, cfW, I, output);
+  Eigen::VectorXi isConstrained = Eigen::VectorXi::Constant(F.rows(),0);
+  Eigen::MatrixXd cfW = Eigen::MatrixXd::Constant(F.rows(),bc.cols(),0);
+
+  for(unsigned i=0; i<b.size();++i)
+  {
+    isConstrained(b(i)) = 1;
+    cfW.row(b(i)) << bc.row(i);
+  }
+  int n = I.rows();
+  igl::GeneralPolyVectorFieldFinder<Eigen::MatrixXd, Eigen::MatrixXi> pvff(V,F,n);
+  pvff.solve(isConstrained, cfW, I, output);
 }
 
 

include/igl/polyvector_field_matchings.cpp

@@ -13,21 +13,36 @@ IGL_INLINE void igl::polyvector_field_matching(
                                                const Eigen::PlainObjectBase<DerivedV>& e,
                                                bool match_with_curl,
                                                Eigen::PlainObjectBase<DerivedM>& mab,
-                                               Eigen::PlainObjectBase<DerivedM>& mba)
+                                               Eigen::PlainObjectBase<DerivedM>& mba,
+                                               bool is_symmetric)
 {
   // make sure the matching preserve ccw order of the vectors across the edge
   // 1) order vectors in a, ccw  e.g. (0,1,2,3)_a not ccw --> (0,3,2,1)_a ccw
   // 2) order vectors in b, ccw  e.g. (0,1,2,3)_b not ccw --> (0,2,1,3)_b ccw
   // 3) the vectors in b that match the ordered vectors in a (in this case  (0,3,2,1)_a ) must be a circular shift of the ccw ordered vectors in b  - so we need to explicitely check only these matchings to find the best ccw one, there are N of them
   int hN = _ua.cols()/3;
-  int N = 2*hN;
-  Eigen::Matrix<typename DerivedS::Scalar,1,Eigen::Dynamic> ua (1,N*3); ua <<_ua, -_ua;
-  Eigen::Matrix<typename DerivedS::Scalar,1,Eigen::Dynamic> ub (1,N*3); ub <<_ub, -_ub;
+  int N;
+  if (is_symmetric)
+    N = 2*hN;
+  else
+    N = hN;
+
+  Eigen::Matrix<typename DerivedS::Scalar,1,Eigen::Dynamic> ua (1,N*3);
+  Eigen::Matrix<typename DerivedS::Scalar,1,Eigen::Dynamic> ub (1,N*3);
+  if (is_symmetric)
+  {
+    ua <<_ua, -_ua;
+    ub <<_ub, -_ub;
+  }
+  else
+  {
+    ua =_ua;
+    ub =_ub;
+  }
+
   Eigen::Matrix<typename DerivedM::Scalar,Eigen::Dynamic,1> order_a, order_b;
-  Eigen::Matrix<typename DerivedS::Scalar, 1, Eigen::Dynamic> sorted_unused;
-  Eigen::VectorXi inv_order_unused;
-  igl::sort_vectors_ccw(ua, na, order_a,false,sorted_unused,false,inv_order_unused);
-  igl::sort_vectors_ccw(ub, nb, order_b,false,sorted_unused,false,inv_order_unused);
+  igl::sort_vectors_ccw(ua, na, order_a);
+  igl::sort_vectors_ccw(ub, nb, order_b);
 
   //checking all possible circshifts of order_b as matches for order_a
   Eigen::Matrix<typename DerivedM::Scalar,Eigen::Dynamic,Eigen::Dynamic> all_matches(N,N);
@@ -71,14 +86,16 @@ IGL_INLINE void igl::polyvector_field_matching(
   for (int i=0; i< N; ++i)
     mba[mab[i]] = i;
 
+  if (is_symmetric)
+  {
   mab = mab.head(hN);
   mba = mba.head(hN);
-
+}
 }
 
 
-template <typename DerivedS, typename DerivedV, typename DerivedF, typename DerivedE, typename DerivedM, typename DerivedC>
-IGL_INLINE typename DerivedC::Scalar igl::polyvector_field_matchings(
+template <typename DerivedS, typename DerivedV, typename DerivedF, typename DerivedE, typename DerivedM>
+IGL_INLINE void igl::polyvector_field_matchings(
                                                                      const Eigen::PlainObjectBase<DerivedS>& sol3D,
                                                                      const Eigen::PlainObjectBase<DerivedV>&V,
                                                                      const Eigen::PlainObjectBase<DerivedF>&F,
@@ -86,9 +103,9 @@ IGL_INLINE typename DerivedC::Scalar igl::polyvector_field_matchings(
                                                                      const Eigen::PlainObjectBase<DerivedV>& FN,
                                                                      const Eigen::MatrixXi &E2F,
                                                                      bool match_with_curl,
+                                                bool is_symmetric,
                                                                      Eigen::PlainObjectBase<DerivedM>& match_ab,
-                                                                     Eigen::PlainObjectBase<DerivedM>& match_ba,
-                                                                     Eigen::PlainObjectBase<DerivedC>& curl)
+                                                Eigen::PlainObjectBase<DerivedM>& match_ba)
 {
   int numEdges = E.rows();
   int half_degree = sol3D.cols()/3;
@@ -101,11 +118,9 @@ IGL_INLINE typename DerivedC::Scalar igl::polyvector_field_matchings(
       isBorderEdge[i] = 1;
   }
 
-  curl.setZero(numEdges,1);
   match_ab.setZero(numEdges, half_degree);
   match_ba.setZero(numEdges, half_degree);
 
-  typename DerivedC::Scalar meanCurl = 0;
   for (int ei=0; ei<numEdges; ++ei)
   {
     if (isBorderEdge[ei])
@@ -124,18 +139,69 @@ IGL_INLINE typename DerivedC::Scalar igl::polyvector_field_matchings(
                                    ce,
                                    match_with_curl,
                                    mab,
-                                   mba);
+                                   mba,
+                                   is_symmetric);
 
     match_ab.row(ei) = mab;
     match_ba.row(ei) = mba;
+  }
+}
+
+
+template <typename DerivedS, typename DerivedV, typename DerivedF, typename DerivedE, typename DerivedM, typename DerivedC>
+IGL_INLINE typename DerivedC::Scalar igl::polyvector_field_matchings(
+                                                                     const Eigen::PlainObjectBase<DerivedS>& sol3D,
+                                                                     const Eigen::PlainObjectBase<DerivedV>&V,
+                                                                     const Eigen::PlainObjectBase<DerivedF>&F,
+                                                                     const Eigen::PlainObjectBase<DerivedE>&E,
+                                                                     const Eigen::PlainObjectBase<DerivedV>& FN,
+                                                                     const Eigen::MatrixXi &E2F,
+                                                                     bool match_with_curl,
+                                                                     bool is_symmetric,
+                                                                     Eigen::PlainObjectBase<DerivedM>& match_ab,
+                                                                     Eigen::PlainObjectBase<DerivedM>& match_ba,
+                                                                     Eigen::PlainObjectBase<DerivedC>& curl)
+{
+  int numEdges = E.rows();
+  int half_degree = sol3D.cols()/3;
+
+  Eigen::VectorXi isBorderEdge;
+  isBorderEdge.setZero(numEdges,1);
+  for(unsigned i=0; i<numEdges; ++i)
+  {
+    if ((E2F(i,0) == -1) || ((E2F(i,1) == -1)))
+      isBorderEdge[i] = 1;
+  }
+
+  igl::polyvector_field_matchings(sol3D, V, F, E, FN, E2F, match_with_curl, is_symmetric, match_ab, match_ba);
+  curl.setZero(numEdges,1);
+  typename DerivedC::Scalar meanCurl = 0;
+  for (int ei=0; ei<numEdges; ++ei)
+  {
+    if (isBorderEdge[ei])
+      continue;
+    // the two faces of the flap
+    int a = E2F(ei,0);
+    int b = E2F(ei,1);
+
+    const Eigen::Matrix<typename DerivedM::Scalar, 1, Eigen::Dynamic> &mab = match_ab.row(ei);
+    const Eigen::Matrix<typename DerivedM::Scalar, 1, Eigen::Dynamic> &mba = match_ba.row(ei);
+
     Eigen::Matrix<typename DerivedS::Scalar, 1, Eigen::Dynamic> matched;
     matched.resize(1, 3*half_degree);
     for (int i = 0; i<half_degree; ++i)
     {
-      int sign = (mab[i]<half_degree)?1:-1;
-      matched.segment(i*3, 3) = sign*sol3D.row(b).segment((mab[i]%half_degree)*3, 3);
+      int sign = 1;
+      int m = mab[i];
+      if (is_symmetric)
+      {
+        sign = (mab[i]<half_degree)?1:-1;
+        m = m%half_degree;
+    }
+      matched.segment(i*3, 3) = sign*sol3D.row(b).segment(m*3, 3);
     }
 
+    Eigen::Matrix<typename DerivedV::Scalar, 1, Eigen::Dynamic> ce = (V.row(E(ei,1))-V.row(E(ei,0))).normalized().template cast<typename DerivedV::Scalar>();
     typename DerivedC::Scalar avgCurl = 0;
     for (int i = 0; i<half_degree; ++i)
       avgCurl += fabs(sol3D.row(a).segment(i*3, 3).dot(ce) - matched.segment(i*3, 3).dot(ce));
@@ -153,14 +219,13 @@ IGL_INLINE typename DerivedC::Scalar igl::polyvector_field_matchings(
   return meanCurl;
 }
 
-
-
 template <typename DerivedS, typename DerivedV, typename DerivedF, typename DerivedM, typename DerivedC>
 IGL_INLINE typename DerivedC::Scalar igl::polyvector_field_matchings(
                                                                      const Eigen::PlainObjectBase<DerivedS>& sol3D,
                                                                      const Eigen::PlainObjectBase<DerivedV>&V,
                                                                      const Eigen::PlainObjectBase<DerivedF>&F,
                                                                      bool match_with_curl,
+                                                                     bool is_symmetric,
                                                                      Eigen::PlainObjectBase<DerivedM>& match_ab,
                                                                      Eigen::PlainObjectBase<DerivedM>& match_ba,
                                                                      Eigen::PlainObjectBase<DerivedC>& curl)
@@ -172,12 +237,31 @@ IGL_INLINE typename DerivedC::Scalar igl::polyvector_field_matchings(
   DerivedV FN;
   igl::per_face_normals(V,F,FN);
 
-  return igl::polyvector_field_matchings(sol3D, V, F, E, FN, E2F, match_with_curl, match_ab, match_ba, curl);
+  return igl::polyvector_field_matchings(sol3D, V, F, E, FN, E2F, match_with_curl, is_symmetric, match_ab, match_ba, curl);
+}
+
+template <typename DerivedS, typename DerivedV, typename DerivedF, typename DerivedM>
+IGL_INLINE void igl::polyvector_field_matchings(
+                                                const Eigen::PlainObjectBase<DerivedS>& sol3D,
+                                                const Eigen::PlainObjectBase<DerivedV>&V,
+                                                const Eigen::PlainObjectBase<DerivedF>&F,
+                                                bool match_with_curl,
+                                                bool is_symmetric,
+                                                Eigen::PlainObjectBase<DerivedM>& match_ab,
+                                                Eigen::PlainObjectBase<DerivedM>& match_ba)
+{
+  Eigen::MatrixXi E, E2F, F2E;
+  igl::edge_topology(V,F,E,F2E,E2F);
+
+  DerivedV FN;
+  igl::per_face_normals(V,F,FN);
+
+  igl::polyvector_field_matchings(sol3D, V, F, E, FN, E2F, match_with_curl, is_symmetric, match_ab, match_ba);
 }
 
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-template Eigen::Matrix<float, -1, 1, 0, -1, 1>::Scalar igl::polyvector_field_matchings<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<float, -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> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 1, 0, -1, 1> >&);
-template Eigen::Matrix<double, -1, 1, 0, -1, 1>::Scalar igl::polyvector_field_matchings<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::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, bool, 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> >&);
+template void igl::polyvector_field_matchings<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::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, bool, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template Eigen::Matrix<double, -1, 1, 0, -1, 1>::Scalar igl::polyvector_field_matchings<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::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, bool, bool, 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> >&);
 #endif

include/igl/polyvector_field_matchings.h

@@ -25,6 +25,10 @@ namespace igl {
   //   e                1 by 3, the vector corresponding to the shared edge between a and b
   //   match_with_curl  boolean flag, determines whether a curl or a smoothness matching will
   //                    be computed
+  //   is_symmetric     boolean flag, determines whether the input vector set field is symmetric(
+  //                    =consists of pairwise collinear vectors in each set, in which case only one
+  //                    of the vectors in the pair is stored) or not, i.e. the set contains all the vectors
+  // )
   // Outputs:
   //   mab              1 by N row vector, describing the matching a->b (i.e. vector #i of the
   //                    vector set in a is matched to vector #mab[i] in b)
@@ -40,7 +44,8 @@ namespace igl {
                                             const Eigen::PlainObjectBase<DerivedV>& e,
                                             bool match_with_curl,
                                             Eigen::PlainObjectBase<DerivedM>& mab,
-                                            Eigen::PlainObjectBase<DerivedM>& mba);
+                                            Eigen::PlainObjectBase<DerivedM>& mba,
+                                            bool is_symmetric);
 
 
   // Given a mesh and a vector set field consisting of unordered N-vector sets defined
@@ -58,6 +63,9 @@ namespace igl {
   //                    via igl::edge_topology)
   //   match_with_curl  boolean flag, determines whether curl or smoothness matchings will
   //                    be computed
+  //   is_symmetric     boolean flag, determines whether the input vector set field is symmetric(
+  //                    =consists of pairwise collinear vectors in each set, in which case only one
+  //                    of the vectors in the pair is stored) or not, i.e. the set contains all the vectors
   // Outputs:
   //   match_ab         #E by N matrix, describing for each edge the matching a->b, where a
   //                    and b are the faces adjacent to the edge (i.e. vector #i of
@@ -78,6 +86,7 @@ namespace igl {
                                                                   const Eigen::PlainObjectBase<DerivedV>& FN,
                                                                   const Eigen::MatrixXi &E2F,
                                                                   bool match_with_curl,
+                                                                  bool is_symmetric,
                                                                   Eigen::PlainObjectBase<DerivedM>& match_ab,
                                                                   Eigen::PlainObjectBase<DerivedM>& match_ba,
                                                                   Eigen::PlainObjectBase<DerivedC>& curl);
@@ -90,10 +99,34 @@ namespace igl {
                                                                   const Eigen::PlainObjectBase<DerivedV>&V,
                                                                   const Eigen::PlainObjectBase<DerivedF>&F,
                                                                   bool match_with_curl,
+                                                                  bool is_symmetric,
                                                                   Eigen::PlainObjectBase<DerivedM>& match_ab,
                                                                   Eigen::PlainObjectBase<DerivedM>& match_ba,
                                                                   Eigen::PlainObjectBase<DerivedC>& curl);
 
+  //Wrappers with no curl output
+  template <typename DerivedS, typename DerivedV, typename DerivedF, typename DerivedM>
+  IGL_INLINE void polyvector_field_matchings(
+                                             const Eigen::PlainObjectBase<DerivedS>& sol3D,
+                                             const Eigen::PlainObjectBase<DerivedV>&V,
+                                             const Eigen::PlainObjectBase<DerivedF>&F,
+                                             bool match_with_curl,
+                                             bool is_symmetric,
+                                             Eigen::PlainObjectBase<DerivedM>& match_ab,
+                                             Eigen::PlainObjectBase<DerivedM>& match_ba);
+  template <typename DerivedS, typename DerivedV, typename DerivedF, typename DerivedE, typename DerivedM>
+  IGL_INLINE void polyvector_field_matchings(
+                                             const Eigen::PlainObjectBase<DerivedS>& sol3D,
+                                             const Eigen::PlainObjectBase<DerivedV>&V,
+                                             const Eigen::PlainObjectBase<DerivedF>&F,
+                                             const Eigen::PlainObjectBase<DerivedE>&E,
+                                             const Eigen::PlainObjectBase<DerivedV>& FN,
+                                             const Eigen::MatrixXi &E2F,
+                                             bool match_with_curl,
+                                             bool is_symmetric,
+                                             Eigen::PlainObjectBase<DerivedM>& match_ab,
+                                             Eigen::PlainObjectBase<DerivedM>& match_ba);
+  
 };
 
 

include/igl/polyvector_field_poisson_reconstruction.cpp

@@ -9,14 +9,12 @@
 
 #include <Eigen/Sparse>
 
-template <typename DerivedV, typename DerivedF, typename DerivedSF, typename DerivedS, typename DerivedE>
-IGL_INLINE double igl::polyvector_field_poisson_reconstruction(
+template <typename DerivedV, typename DerivedF, typename DerivedSF, typename DerivedS>
+IGL_INLINE void igl::polyvector_field_poisson_reconstruction(
   const Eigen::PlainObjectBase<DerivedV> &Vcut,
   const Eigen::PlainObjectBase<DerivedF> &Fcut,
   const Eigen::PlainObjectBase<DerivedS> &sol3D_combed,
-  Eigen::PlainObjectBase<DerivedSF> &scalars,
-  Eigen::PlainObjectBase<DerivedS> &sol3D_recon,
-  Eigen::PlainObjectBase<DerivedE> &max_error )
+                                                               Eigen::PlainObjectBase<DerivedSF> &scalars)
   {
     Eigen::SparseMatrix<typename DerivedV::Scalar> gradMatrix;
     igl::grad(Vcut, Fcut, gradMatrix);
@@ -34,8 +32,6 @@ IGL_INLINE double igl::polyvector_field_poisson_reconstruction(
 
     int half_degree = sol3D_combed.cols()/3;
 
-    sol3D_recon.setZero(sol3D_combed.rows(),sol3D_combed.cols());
-
     int numF = Fcut.rows();
     scalars.setZero(Vcut.rows(),half_degree);
 
@@ -65,13 +61,34 @@ IGL_INLINE double igl::polyvector_field_poisson_reconstruction(
       Eigen::VectorXd bu = igl::slice(b, Iu);
 
       Eigen::VectorXd rhs = bu-Quk*xk;
-      Eigen::MatrixXd yu = solver.solve(rhs);
+      Eigen::VectorXd yu = solver.solve(rhs);
 
-      Eigen::VectorXi index = i*Eigen::VectorXi::Ones(Iu.rows(),1);
-      igl::slice_into(yu, Iu, index, scalars);scalars(Ik[0],i)=xk[0];
+      Eigen::VectorXd y(Vcut.rows(),1);
+      igl::slice_into(yu, Iu, 1, y);y(Ik[0])=xk[0];
+      scalars.col(i) = y;
     }
+}
+
+template <typename DerivedV, typename DerivedF, typename DerivedSF, typename DerivedS, typename DerivedE>
+IGL_INLINE double igl::polyvector_field_poisson_reconstruction(
+                                                               const Eigen::PlainObjectBase<DerivedV> &Vcut,
+                                                               const Eigen::PlainObjectBase<DerivedF> &Fcut,
+                                                               const Eigen::PlainObjectBase<DerivedS> &sol3D_combed,
+                                                               Eigen::PlainObjectBase<DerivedSF> &scalars,
+                                                               Eigen::PlainObjectBase<DerivedS> &sol3D_recon,
+                                                               Eigen::PlainObjectBase<DerivedE> &max_error )
+{
+  
+  igl::polyvector_field_poisson_reconstruction(Vcut, Fcut, sol3D_combed, scalars);
+
+  Eigen::SparseMatrix<typename DerivedV::Scalar> gradMatrix;
+  igl::grad(Vcut, Fcut, gradMatrix);
+  int numF = Fcut.rows();
+  int half_degree = sol3D_combed.cols()/3;
 
     //    evaluate gradient of found scalar function
+  sol3D_recon.setZero(sol3D_combed.rows(),sol3D_combed.cols());
+  
     for (int i =0; i<half_degree; ++i)
     {
       Eigen::VectorXd vec_poisson = gradMatrix*scalars.col(i);
@@ -91,7 +108,9 @@ IGL_INLINE double igl::polyvector_field_poisson_reconstruction(
     return max_error.mean();
   }
 
+
   #ifdef IGL_STATIC_LIBRARY
   // Explicit template specialization
   template double igl::polyvector_field_poisson_reconstruction<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::Matrix<float, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -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<float, -1, 1, 0, -1, 1> >&);
+template double igl::polyvector_field_poisson_reconstruction<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::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -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<double, -1, 1, 0, -1, 1> >&);
   #endif

include/igl/polyvector_field_poisson_reconstruction.h

@@ -41,6 +41,13 @@ namespace igl {
                                                             Eigen::PlainObjectBase<DerivedSF> &scalars,
                                                             Eigen::PlainObjectBase<DerivedS> &sol3D_recon,
                                                             Eigen::PlainObjectBase<DerivedE> &max_error );
+  //Wrappers with less output
+  template <typename DerivedV, typename DerivedF, typename DerivedSF, typename DerivedS>
+  IGL_INLINE void polyvector_field_poisson_reconstruction(
+                                                            const Eigen::PlainObjectBase<DerivedV> &Vcut,
+                                                            const Eigen::PlainObjectBase<DerivedF> &Fcut,
+                                                            const Eigen::PlainObjectBase<DerivedS> &sol3D_combed,
+                                                            Eigen::PlainObjectBase<DerivedSF> &scalars);
 
 
 };

include/igl/polyvector_field_singularities_from_matchings.cpp

@@ -1,3 +1,4 @@
+#include <iostream>
 #include <igl/polyvector_field_singularities_from_matchings.h>
 #include <igl/is_border_vertex.h>
 #include <igl/vertex_triangle_adjacency.h>
@@ -15,52 +16,86 @@ void igl::polyvector_field_one_ring_matchings(const Eigen::PlainObjectBase<Deriv
                                               const Eigen::PlainObjectBase<DerivedM> &match_ab,
                                               const Eigen::PlainObjectBase<DerivedM> &match_ba,
                                               const int vi,
-                                              const int vector_to_match,
-                                              Eigen::VectorXi &mvi,
+                                              Eigen::MatrixXi &mvi,
                                               Eigen::VectorXi &fi)
 {
   int half_degree = match_ab.cols();
-  mvi.resize(VF[vi].size()+1,1);
+  mvi.resize(VF[vi].size()+1,half_degree);
   fi.resize(VF[vi].size()+1,1);
   //start from one face
-  const int &fstart = VF[vi][0];
+  //first, check if the vertex is on a boundary
+  //then there must be two faces that are on the boundary
+  //(other cases not supported)
+
+  int fstart = -1;
+  int ind = 0;
+  for (int i =0; i<VF[vi].size(); ++i)
+  {
+    int fi = VF[vi][i];
+    for (int  j=0; j<3; ++j)
+      if (F(fi,j)==vi && TT(fi,j) == -1)
+      {
+        ind ++;
+        fstart = fi;
+        //        break;
+      }
+  }
+  if (ind >1 )
+  {
+    std::cerr<<"igl::polyvector_field_one_ring_matchings -- vertex "<<vi<< " is on an unusual boundary"<<std::endl;
+    exit(1);
+  }
+  if (fstart == -1)
+    fstart = VF[vi][0];
   int current_face = fstart;
   int i =0;
-  mvi[i] = vector_to_match;
   fi[i] = current_face;
+  for (int j=0; j<half_degree; ++j)
+    mvi(i,j) = j;
+
   int next_face = -1;
-  while (next_face != fstart)
+  while (next_face != fstart && current_face!=-1)
   {
     // look for the vertex
     int j=-1;
     for (unsigned z=0; z<3; ++z)
-      if (F(current_face,z) == vi)
+      if (F(current_face,(z+1)%3) == vi)
+      {
         j=z;
+        break;
+      }
     assert(j!=-1);
 
     next_face = TT(current_face, j);
     ++i;
 
-    // look at the edge between the two faces
-    const int &current_edge = F2E(current_face,j);
-
-    // check its orientation to determine whether match_ab or match_ba should be used
-    if ((E2F(current_edge,0) == current_face) &&
-        (E2F(current_edge,1) == next_face) )
-    {
-      //look at match_ab
-      mvi[i] = match_ab(current_edge,(mvi[i-1])%half_degree);
-    }
+    if (next_face == -1)
+      mvi.row(i).setConstant(-1);
     else
     {
-      assert((E2F(current_edge,1) == current_face) &&
-             (E2F(current_edge,0) == next_face));
-      //look at match_ba
-      mvi[i] = match_ba(current_edge,(mvi[i-1])%half_degree);
+      // look at the edge between the two faces
+      const int &current_edge = F2E(current_face,j);
+
+      for (int k=0; k<half_degree; ++k)
+      {
+        // check its orientation to determine whether match_ab or match_ba should be used
+        if ((E2F(current_edge,0) == current_face) &&
+            (E2F(current_edge,1) == next_face) )
+        {
+          //look at match_ab
+          mvi(i,k) = match_ab(current_edge,(mvi(i-1,k))%half_degree);
+        }
+        else
+        {
+          assert((E2F(current_edge,1) == current_face) &&
+                 (E2F(current_edge,0) == next_face));
+          //look at match_ba
+          mvi(i,k) = match_ba(current_edge,(mvi(i-1,k))%half_degree);
+        }
+        if (mvi(i-1,k)>=half_degree)
+          mvi(i,k) = (mvi(i,k)+half_degree)%(2*half_degree);
+      }
     }
-    if (mvi[i-1]>=half_degree)
-      mvi[i] = (mvi[i]+half_degree)%(2*half_degree);
-
     current_face = next_face;
     fi[i] = current_face;
   }
@@ -106,20 +141,24 @@ IGL_INLINE void igl::polyvector_field_singularities_from_matchings(
 
   std::vector<int> singularities_v;
   int half_degree = match_ab.cols();
-  //pick one of the vectors to check for singularities
-  for (int vector_to_match = 0; vector_to_match < half_degree; ++vector_to_match)
+  for (int vi =0; vi<numV; ++vi)
   {
-    for (int vi =0; vi<numV; ++vi)
+    ///check that is on border..
+    if (V_border[vi])
+      continue;
+    Eigen::VectorXi fi;
+    Eigen::MatrixXi mvi;
+    igl::polyvector_field_one_ring_matchings(V, F, VF, E2F, F2E, TT, match_ab, match_ba, vi, mvi, fi);
+
+    int num = fi.size();
+    //pick one of the vectors to check for singularities
+    for (int vector_to_match = 0; vector_to_match < half_degree; ++vector_to_match)
     {
-      ///check that is on border..
-      if (V_border[vi])
-        continue;
-      Eigen::VectorXi mvi,fi;
-
-      igl::polyvector_field_one_ring_matchings(V, F, VF, E2F, F2E, TT, match_ab, match_ba, vi, vector_to_match, mvi, fi);
-
-      if(mvi.tail(1) != mvi.head(1))
+      if(mvi(num-1,vector_to_match) != mvi(0,vector_to_match))
+      {
         singularities_v.push_back(vi);
+        break;
+      }
     }
   }
   std::sort(singularities_v.begin(), singularities_v.end());
@@ -129,8 +168,66 @@ IGL_INLINE void igl::polyvector_field_singularities_from_matchings(
   igl::list_to_matrix(singularities_v, singularities);
 }
 
+
+template <typename DerivedV, typename DerivedF, typename DerivedM, typename DerivedS>
+IGL_INLINE void igl::polyvector_field_singularities_from_matchings(
+                                                                   const Eigen::PlainObjectBase<DerivedV> &V,
+                                                                   const Eigen::PlainObjectBase<DerivedF> &F,
+                                                                   const Eigen::PlainObjectBase<DerivedM> &match_ab,
+                                                                   const Eigen::PlainObjectBase<DerivedM> &match_ba,
+                                                                   Eigen::PlainObjectBase<DerivedS> &singularities,
+                                                                   Eigen::PlainObjectBase<DerivedS> &singularity_indices)
+{
+
+  std::vector<bool> V_border = igl::is_border_vertex(V,F);
+  std::vector<std::vector<int> > VF, VFi;
+  igl::vertex_triangle_adjacency(V,F,VF,VFi);
+
+  Eigen::MatrixXi TT, TTi;
+  igl::triangle_triangle_adjacency(V,F,TT,TTi);
+
+  Eigen::MatrixXi E, E2F, F2E;
+  igl::edge_topology(V,F,E,F2E,E2F);
+
+  igl::polyvector_field_singularities_from_matchings(V, F, V_border, VF, TT, E2F, F2E, match_ab, match_ba, singularities, singularity_indices);
+}
+
+template <typename DerivedV, typename DerivedF, typename DerivedM, typename VFType, typename DerivedS>
+IGL_INLINE void igl::polyvector_field_singularities_from_matchings(
+                                                                   const Eigen::PlainObjectBase<DerivedV> &V,
+                                                                   const Eigen::PlainObjectBase<DerivedF> &F,
+                                                                   const std::vector<bool> &V_border,
+                                                                   const std::vector<std::vector<VFType> > &VF,
+                                                                   const Eigen::MatrixXi &TT,
+                                                                   const Eigen::MatrixXi &E2F,
+                                                                   const Eigen::MatrixXi &F2E,
+                                                                   const Eigen::PlainObjectBase<DerivedM> &match_ab,
+                                                                   const Eigen::PlainObjectBase<DerivedM> &match_ba,
+                                                                   Eigen::PlainObjectBase<DerivedS> &singularities,
+                                                                   Eigen::PlainObjectBase<DerivedS> &singularity_indices)
+{
+  igl::polyvector_field_singularities_from_matchings(V, F, V_border, VF, TT, E2F, F2E, match_ab, match_ba, singularities);
+
+  singularity_indices.setZero(singularities.size(), 1);
+
+  //get index from first vector only
+  int vector_to_match = 0;
+  for (int i =0; i<singularities.size(); ++i)
+  {
+    int vi = singularities[i];
+
+    // Eigen::VectorXi mvi,fi;
+    // igl::polyvector_field_one_ring_matchings(V, F, VF, E2F, F2E, TT, match_ab, match_ba, vi, vector_to_match, mvi, fi);
+    Eigen::VectorXi fi;
+    Eigen::MatrixXi mvi;
+    igl::polyvector_field_one_ring_matchings(V, F, VF, E2F, F2E, TT, match_ab, match_ba, vi, mvi, fi);
+
+    singularity_indices[i] = (mvi(mvi.rows()-1,vector_to_match) - vector_to_match);
+  }
+
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-template void igl::polyvector_field_singularities_from_matchings<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, int, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<bool, std::allocator<bool> > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::polyvector_field_singularities_from_matchings<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<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/polyvector_field_singularities_from_matchings.h

@@ -17,7 +17,7 @@ namespace igl {
   
   // We are given a polyvector field on a mesh (with N vectors per face), and matchings between the vector sets
   // across all non-boundary edges.  The function computes, for the one-ring
-  // neighborhood of a given vertex, and for the selected vector of the vector set,
+  // neighborhood of a given vertex, and for each of the vectors of the vector set,
   // the sequence of the vectors that match to it around the one-ring. If the vector that
   // we land on by following the matchings is not the original vector that we started from,
   // the vertex is a singularity.
@@ -40,13 +40,13 @@ namespace igl {
   //                    and b are the faces adjacent to the edge (i.e. vector #mba[i] of
   //                    the vector set in a is matched to vector #i in b)
   //   vi               the selected one ring
-  //   vector_id        the selected vector of the polyvector
   //
   // Output:
   //   mvi              #numOneRingFaces by 1 list of the indices of the sequentially matching vectors
   //                    in the faces of the one ring (first enty is always vector_id, then the vector matching
   //                    vector_id in the next face, then the vector matching that in the third face etc.)
-  //   fi               #numOneRingFaces by 1 list of the sequentially visited faces in the one ring neighborhood
+  //   fi               #numOneRingFaces by 1 list of the sequentially visited faces in the one ring neighborhood.
+  //                    The one-ring is traversed in CLOCKWISE order with respect to the outward normal. (=opposite)
   //
   template <typename DerivedV, typename DerivedF, typename DerivedM, typename VFType, typename DerivedTT>
   void polyvector_field_one_ring_matchings(const Eigen::PlainObjectBase<DerivedV> &V,
@@ -58,8 +58,7 @@ namespace igl {
                                            const Eigen::PlainObjectBase<DerivedM> &match_ab,
                                            const Eigen::PlainObjectBase<DerivedM> &match_ba,
                                            const int vi,
-                                           const int vector_id,
-                                           Eigen::VectorXi &mvi,
+                                           Eigen::MatrixXi &mvi,
                                            Eigen::VectorXi &fi);
   
   
@@ -114,6 +113,31 @@ namespace igl {
                                                                 Eigen::PlainObjectBase<DerivedS> &singularities);
   
   
+  // Same pair as above but also returns singularity indices
+  template <typename DerivedV, typename DerivedF, typename DerivedM, typename DerivedS>
+  IGL_INLINE void polyvector_field_singularities_from_matchings(
+                                                                const Eigen::PlainObjectBase<DerivedV> &V,
+                                                                const Eigen::PlainObjectBase<DerivedF> &F,
+                                                                const Eigen::PlainObjectBase<DerivedM> &match_ab,
+                                                                const Eigen::PlainObjectBase<DerivedM> &match_ba,
+                                                                Eigen::PlainObjectBase<DerivedS> &singularities,
+                                                                Eigen::PlainObjectBase<DerivedS> &singularity_indices);
+  
+  template <typename DerivedV, typename DerivedF, typename DerivedM, typename VFType, typename DerivedS>
+  IGL_INLINE void polyvector_field_singularities_from_matchings(
+                                                                const Eigen::PlainObjectBase<DerivedV> &V,
+                                                                const Eigen::PlainObjectBase<DerivedF> &F,
+                                                                const std::vector<bool> &V_border,
+                                                                const std::vector<std::vector<VFType> > &VF,
+                                                                const Eigen::MatrixXi &TT,
+                                                                const Eigen::MatrixXi &E2F,
+                                                                const Eigen::MatrixXi &F2E,
+                                                                const Eigen::PlainObjectBase<DerivedM> &match_ab,
+                                                                const Eigen::PlainObjectBase<DerivedM> &match_ba,
+                                                                Eigen::PlainObjectBase<DerivedS> &singularities,
+                                                                Eigen::PlainObjectBase<DerivedS> &singularity_indices);
+
+  
   
 };
 

include/igl/rotate_vectors.cpp

@@ -17,6 +17,8 @@ IGL_INLINE Eigen::MatrixXd igl::rotate_vectors(
 
   for (unsigned i=0; i<V.rows();++i)
   {
+    double norm = V.row(i).norm();
+    
     // project onto the tangent plane and convert to angle
     double a = atan2(B2.row(i).dot(V.row(i)),B1.row(i).dot(V.row(i)));
 
@@ -24,7 +26,7 @@ IGL_INLINE Eigen::MatrixXd igl::rotate_vectors(
     a += (A.size() == 1) ? A(0) : A(i);
 
     // move it back to global coordinates
-    RV.row(i) = cos(a) * B1.row(i) + sin(a) * B2.row(i);
+    RV.row(i) = norm*cos(a) * B1.row(i) + norm*sin(a) * B2.row(i);
   }
 
   return RV;

include/igl/seam_edges.cpp

@@ -0,0 +1,180 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Yotam Gingold <yotam@yotamgingold.com>
+// 
+// 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 "seam_edges.h"
+#include <unordered_map>
+#include <unordered_set>
+#include <cassert>
+
+namespace {
+    // Computes the orientation of `c` relative to the line between `a` and `b`.
+    // Assumes 2D vector input.
+    // Based on: https://www.cs.cmu.edu/~quake/robust.html
+    template< typename scalar_t >
+    inline scalar_t orientation(
+        const Eigen::Matrix< scalar_t, 2, 1 >& a,
+        const Eigen::Matrix< scalar_t, 2, 1 >& b,
+        const Eigen::Matrix< scalar_t, 2, 1 >& c
+        ) {
+        typedef Eigen::Matrix< scalar_t, 2, 1 > Vector2S;
+        const Vector2S row0 = a - c;
+        const Vector2S row1 = b - c;
+        return row0(0)*row1(1) - row1(0)*row0(1);
+    }
+}
+
+// I have verified that this function produces the exact same output as
+// `find_seam_fast.py` for `cow_triangled.obj`.
+template <typename DerivedV, typename DerivedF, typename DerivedT>
+IGL_INLINE void seam_edges(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedT>& TC,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    const Eigen::PlainObjectBase<DerivedF>& FTC,
+    Eigen::PlainObjectBase<DerivedF>& seams,
+    Eigen::PlainObjectBase<DerivedF>& boundaries,
+    Eigen::PlainObjectBase<DerivedF>& foldovers
+    )
+{
+    // Assume triangles.
+    assert( F.cols() == 3 );
+    assert( F.cols() == FTC.cols() );
+    assert( F.rows() == FTC.rows() );
+    
+    // Assume 2D texture coordinates (foldovers tests).
+    assert( TC.cols() == 2 );
+    typedef Eigen::Matrix< typename DerivedT::Scalar, 2, 1 > Vector2S;
+    
+    seams     .setZero( 3*F.rows(), 4 );
+    boundaries.setZero( 3*F.rows(), 2 );
+    foldovers .setZero( 3*F.rows(), 4 );
+    
+    int num_seams = 0;
+    int num_boundaries = 0;
+    int num_foldovers = 0;
+    
+    // A map from a pair of vertex indices to the index (face and endpoints)
+    // into face_position_indices.
+    // The following should be true for every key, value pair:
+    //    key == face_position_indices[ value ]
+    // This gives us a "reverse map" so that we can look up other face
+    // attributes based on position edges.
+    // The value are written in the format returned by numpy.where(),
+    // which stores multi-dimensional indices such as array[a0,b0], array[a1,b1]
+    // as ( (a0,a1), (b0,b1) ).
+    
+    // We need to make a hash function for our directed edges.
+    // We'll use i*V.rows() + j.
+    typedef std::pair< typename DerivedF::Scalar, typename DerivedF::Scalar > directed_edge;
+	const int numV = V.rows();
+	const int numF = F.rows();
+	auto edge_hasher = [numV]( directed_edge const& e ) { return e.first*numV + e.second; };
+	// When we pass a hash function object, we also need to specify the number of buckets.
+	// The Euler characteristic says that the number of undirected edges is numV + numF -2*genus.
+	std::unordered_map< directed_edge, std::pair< int, int >, decltype( edge_hasher ) > directed_position_edge2face_position_index( 2*( numV + numF ), edge_hasher );
+    for( int fi = 0; fi < F.rows(); ++fi ) {
+        for( int i = 0; i < 3; ++i ) {
+            const int j = ( i+1 ) % 3;
+            directed_position_edge2face_position_index[ std::make_pair( F(fi,i), F(fi,j) ) ] = std::make_pair( fi, i );
+        }
+    }
+    
+    // First find all undirected position edges (collect both orientations of
+    // the directed edges).
+    std::unordered_set< directed_edge, decltype( edge_hasher ) > undirected_position_edges( numV + numF, edge_hasher );
+    for( const auto& el : directed_position_edge2face_position_index ) {
+        undirected_position_edges.insert( el.first );
+        undirected_position_edges.insert( std::make_pair( el.first.second, el.first.first ) );
+    }
+    
+    // Now we will iterate over all position edges.
+    // Seam edges are the edges whose two opposite directed edges have different
+    // texcoord indices (or one doesn't exist at all in the case of a mesh
+    // boundary).
+    for( const auto& vp_edge : undirected_position_edges ) {
+        const auto vp_edge_reverse = std::make_pair( vp_edge.second, vp_edge.first );
+        // If it and its opposite exist as directed edges, check if their
+        // texture coordinate indices match.
+        if( directed_position_edge2face_position_index.count( vp_edge ) &&
+            directed_position_edge2face_position_index.count( vp_edge_reverse ) ) {
+            const auto forwards = directed_position_edge2face_position_index[ vp_edge ];
+            const auto backwards = directed_position_edge2face_position_index[ vp_edge_reverse ];
+            
+            // NOTE: They should never be equal.
+            assert( forwards != backwards );
+            
+            // NOTE: Non-matching seam edges will show up twice, once as
+            //       ( forwards, backwards ) and once as
+            //       ( backwards, forwards ). We don't need to examine both,
+            //       so continue only if forwards < backwards.
+            if( forwards < backwards ) continue;
+
+            // If the texcoord indices match (are similarly flipped),
+            // this edge is not a seam. It could be a foldover.
+            if( std::make_pair( FTC( forwards.first, forwards.second ), FTC( forwards.first, ( forwards.second+1 ) % 3 ) )
+                ==
+                std::make_pair( FTC( backwards.first, ( backwards.second+1 ) % 3 ), FTC( backwards.first, backwards.second ) )
+                ) {
+                
+                // Check for foldovers in UV space.
+                // Get the edge (a,b) and the two opposite vertices's texture coordinates.
+                const Vector2S a = TC.row( FTC( forwards.first,  forwards.second ) );
+                const Vector2S b = TC.row( FTC( forwards.first, (forwards.second+1) % 3 ) );
+                const Vector2S c_forwards  = TC.row( FTC( forwards .first, (forwards .second+2) % 3 ) );
+                const Vector2S c_backwards = TC.row( FTC( backwards.first, (backwards.second+2) % 3 ) );
+                // If the opposite vertices' texture coordinates fall on the same side
+                // of the edge, we have a UV-space foldover.
+                const auto orientation_forwards = orientation( a, b, c_forwards );
+                const auto orientation_backwards = orientation( a, b, c_backwards );
+                if( ( orientation_forwards > 0 && orientation_backwards > 0 ) ||
+                    ( orientation_forwards < 0 && orientation_backwards < 0 )
+                    ) {
+                    foldovers( num_foldovers, 0 ) = forwards.first;
+                    foldovers( num_foldovers, 1 ) = forwards.second;
+                    foldovers( num_foldovers, 2 ) = backwards.first;
+                    foldovers( num_foldovers, 3 ) = backwards.second;
+                    num_foldovers += 1;
+                }
+            }
+            
+            // Otherwise, we have a non-matching seam edge.
+            else {
+                seams( num_seams, 0 ) = forwards.first;
+                seams( num_seams, 1 ) = forwards.second;
+                seams( num_seams, 2 ) = backwards.first;
+                seams( num_seams, 3 ) = backwards.second;
+                num_seams += 1;
+            }
+        }
+        // Otherwise, the edge and its opposite aren't both in the directed
+        // edges. One of them should be.
+        else if( directed_position_edge2face_position_index.count( vp_edge ) ) {
+            const auto forwards = directed_position_edge2face_position_index[ vp_edge ];
+            boundaries( num_boundaries, 0 ) = forwards.first;
+            boundaries( num_boundaries, 1 ) = forwards.second;
+            num_boundaries += 1;
+        }
+        else if( directed_position_edge2face_position_index.count( vp_edge_reverse ) ) {
+            const auto backwards = directed_position_edge2face_position_index[ vp_edge_reverse ];
+            boundaries( num_boundaries, 0 ) = backwards.first;
+            boundaries( num_boundaries, 1 ) = backwards.second;
+            num_boundaries += 1;
+        }
+        else {
+            // This should never happen! One of these two must have been seen.
+            assert(
+                directed_position_edge2face_position_index.count( vp_edge ) ||
+                directed_position_edge2face_position_index.count( vp_edge_reverse )
+                );
+        }
+    }
+    
+    seams     .conservativeResize( num_seams,      Eigen::NoChange_t() );
+    boundaries.conservativeResize( num_boundaries, Eigen::NoChange_t() );
+    foldovers .conservativeResize( num_foldovers,  Eigen::NoChange_t() );
+}

include/igl/seam_edges.h

@@ -0,0 +1,57 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Yotam Gingold <yotam@yotamgingold.com>
+// 
+// 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/.
+#ifndef IGL_SEAM_EDGES_H
+#define IGL_SEAM_EDGES_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+    /*
+    Returns all UV-space boundaries of a mesh.
+    Inputs:
+        V: The positions of the input mesh.
+        TC: The 2D texture coordinates of the input mesh (2 columns).
+        F: A manifold-with-boundary triangle mesh, as vertex indices into `V` for the three vertices of each triangle.
+        FTC: Indices into `TC` for the three vertices of each triangle.
+    Outputs:
+        seams: Edges where the forwards and backwards directions have different texture
+               coordinates, as a #seams-by-4 matrix of indices.
+               Each row is organized as [ forward_face_index, forward_face_vertex_index,
+               backwards_face_index, backwards_face_vertex_index ]
+               such that one side of the seam is the edge:
+                    F[ seams( i, 0 ), seams( i, 1 ) ], F[ seams( i, 0 ), (seams( i, 1 ) + 1) % 3 ]
+               and the other side is the edge:
+                    F[ seams( i, 2 ), seams( i, 3 ) ], F[ seams( i, 2 ), (seams( i, 3 ) + 1) % 3 ]
+        boundaries: Edges with only one incident triangle, as a #boundaries-by-2 matrix of
+                    indices. Each row is organized as [ face_index, face_vertex_index ]
+                    such that the edge is:
+                        F[ boundaries( i, 0 ), boundaries( i, 1 ) ], F[ boundaries( i, 0 ), (boundaries( i, 1 ) + 1) % 3 ]
+        foldovers: Edges where the two incident triangles fold over each other in UV-space,
+                   as a #foldovers-by-4 matrix of indices.
+                   Each row is organized as [ forward_face_index, forward_face_vertex_index,
+                   backwards_face_index, backwards_face_vertex_index ]
+                   such that one side of the foldover is the edge:
+                        F[ foldovers( i, 0 ), foldovers( i, 1 ) ], F[ foldovers( i, 0 ), (foldovers( i, 1 ) + 1) % 3 ]
+                   and the other side is the edge:
+                        F[ foldovers( i, 2 ), foldovers( i, 3 ) ], F[ foldovers( i, 2 ), (foldovers( i, 3 ) + 1) % 3 ]
+    */
+    template <typename DerivedV, typename DerivedF, typename DerivedT>
+    IGL_INLINE void seam_edges(
+        const Eigen::PlainObjectBase<DerivedV>& V,
+        const Eigen::PlainObjectBase<DerivedT>& TC,
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        const Eigen::PlainObjectBase<DerivedF>& FTC,
+        Eigen::PlainObjectBase<DerivedF>& seams,
+        Eigen::PlainObjectBase<DerivedF>& boundaries,
+        Eigen::PlainObjectBase<DerivedF>& foldovers
+        );
+}
+#ifndef IGL_STATIC_LIBRARY
+#   include "seam_edges.cpp"
+#endif
+#endif // IGL_SEAM_EDGES_H

include/igl/slice.cpp

@@ -278,6 +278,9 @@ template void igl::slice<Eigen::SparseMatrix<int, 0, int>, Eigen::Matrix<int, -1
 template Eigen::Matrix<double, -1, -1, 0, -1, -1> igl::slice<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&);
 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> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::slice<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::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::slice<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> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, Eigen::Matrix<double, -1, 1, 0, -1, 1>&);
+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::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, Eigen::Matrix<int, -1, -1, 0, -1, -1>&);
+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::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, Eigen::Matrix<int, -1, 1, 0, -1, 1>&);
 template void igl::slice<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::slice<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> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
 template void igl::slice<Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >, 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> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);

include/igl/sort_vectors_ccw.cpp

@@ -6,14 +6,9 @@ template <typename DerivedS, typename DerivedI>
 IGL_INLINE void igl::sort_vectors_ccw(
   const Eigen::PlainObjectBase<DerivedS>& P,
   const Eigen::PlainObjectBase<DerivedS>& N,
-  Eigen::PlainObjectBase<DerivedI> &order,
-  const bool do_sorted,
-  Eigen::PlainObjectBase<DerivedS> &sorted,
-  const bool do_inv_order,
-  Eigen::PlainObjectBase<DerivedI> &inv_order)
+  Eigen::PlainObjectBase<DerivedI> &order)
 {
   int half_degree = P.cols()/3;
-
   //local frame
   Eigen::Matrix<typename DerivedS::Scalar,1,3> e1 = P.head(3).normalized();
   Eigen::Matrix<typename DerivedS::Scalar,1,3> e3 = N.normalized();
@@ -25,7 +20,7 @@ IGL_INLINE void igl::sort_vectors_ccw(
   for (int i=0; i<half_degree; ++i)
   {
     Eigen::Matrix<typename DerivedS::Scalar,1,3> Pl = F.colPivHouseholderQr().solve(P.segment(i*3,3).transpose()).transpose();
-    assert(fabs(Pl(2))/Pl.cwiseAbs().maxCoeff() <1e-5);
+//    assert(fabs(Pl(2))/Pl.cwiseAbs().maxCoeff() <1e-5);
     angles[i] = atan2(Pl(1),Pl(0));
   }
 
@@ -39,14 +34,31 @@ IGL_INLINE void igl::sort_vectors_ccw(
       order[i] = order[i+1];
     order(half_degree-1) = temp;
   }
-  if (do_sorted)
+}
+
+template <typename DerivedS, typename DerivedI>
+IGL_INLINE void igl::sort_vectors_ccw(
+  const Eigen::PlainObjectBase<DerivedS>& P,
+  const Eigen::PlainObjectBase<DerivedS>& N,
+  Eigen::PlainObjectBase<DerivedI> &order,
+  Eigen::PlainObjectBase<DerivedS> &sorted)
   {
+  int half_degree = P.cols()/3;
+  igl::sort_vectors_ccw(P,N,order);
     sorted.resize(1,half_degree*3);
     for (int i=0; i<half_degree; ++i)
       sorted.segment(i*3,3) = P.segment(order[i]*3,3);
   }
-  if (do_inv_order)
+
+template <typename DerivedS, typename DerivedI>
+IGL_INLINE void igl::sort_vectors_ccw(
+  const Eigen::PlainObjectBase<DerivedS>& P,
+  const Eigen::PlainObjectBase<DerivedS>& N,
+  Eigen::PlainObjectBase<DerivedI> &order,
+  Eigen::PlainObjectBase<DerivedI> &inv_order)
   {
+  int half_degree = P.cols()/3;
+  igl::sort_vectors_ccw(P,N,order);
     inv_order.resize(half_degree,1);
     for (int i=0; i<half_degree; ++i)
     {
@@ -60,9 +72,25 @@ IGL_INLINE void igl::sort_vectors_ccw(
     assert(inv_order[0] == 0);
   }
 
+template <typename DerivedS, typename DerivedI>
+IGL_INLINE void igl::sort_vectors_ccw(
+  const Eigen::PlainObjectBase<DerivedS>& P,
+  const Eigen::PlainObjectBase<DerivedS>& N,
+  Eigen::PlainObjectBase<DerivedI> &order,
+  Eigen::PlainObjectBase<DerivedS> &sorted,
+  Eigen::PlainObjectBase<DerivedI> &inv_order)
+{
+  int half_degree = P.cols()/3;
+
+  igl::sort_vectors_ccw(P,N,order,inv_order);
+
+  sorted.resize(1,half_degree*3);
+  for (int i=0; i<half_degree; ++i)
+    sorted.segment(i*3,3) = P.segment(order[i]*3,3);
 }
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-template void igl::sort_vectors_ccw<Eigen::Matrix<double, 1, -1, 1, 1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> >&, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::sort_vectors_ccw<Eigen::Matrix<double, 1, -1, 1, 1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::sort_vectors_ccw<Eigen::Matrix<double, 1, -1, 1, 1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> >&);
 #endif

include/igl/sort_vectors_ccw.h

@@ -21,8 +21,6 @@ namespace igl {
   // Inputs:
   //   P               1 by 3N row vector of the vectors to be sorted, stacked horizontally
   //   N               #1 by 3 normal of the plane where the vectors lie
-  //   do_sorted       boolean flag, determines whether to return the sorted vector set
-  //   do_inv_order    boolean flag, determines whether to return the "inverse" set of indices
   // Output:
   //   order           N by 1 order of the vectors (indices of the unordered vectors into
   //                   the ordered vector set)
@@ -35,10 +33,30 @@ namespace igl {
                                    const Eigen::PlainObjectBase<DerivedS>& P,
                                    const Eigen::PlainObjectBase<DerivedS>& N,
                                    Eigen::PlainObjectBase<DerivedI> &order,
-                                   const bool do_sorted,
                                    Eigen::PlainObjectBase<DerivedS> &sorted,
-                                   const bool do_inv_order,
                                    Eigen::PlainObjectBase<DerivedI> &inv_order);
+
+   template <typename DerivedS, typename DerivedI>
+   IGL_INLINE void sort_vectors_ccw(
+                                    const Eigen::PlainObjectBase<DerivedS>& P,
+                                    const Eigen::PlainObjectBase<DerivedS>& N,
+                                    Eigen::PlainObjectBase<DerivedI> &order,
+                                    Eigen::PlainObjectBase<DerivedS> &sorted);
+
+    template <typename DerivedS, typename DerivedI>
+    IGL_INLINE void sort_vectors_ccw(
+                                     const Eigen::PlainObjectBase<DerivedS>& P,
+                                     const Eigen::PlainObjectBase<DerivedS>& N,
+                                     Eigen::PlainObjectBase<DerivedI> &order,
+                                     Eigen::PlainObjectBase<DerivedI> &inv_order);
+
+
+     template <typename DerivedS, typename DerivedI>
+     IGL_INLINE void sort_vectors_ccw(
+                                      const Eigen::PlainObjectBase<DerivedS>& P,
+                                      const Eigen::PlainObjectBase<DerivedS>& N,
+                                      Eigen::PlainObjectBase<DerivedI> &order);
+
 };
 
 

include/igl/streamlines.cpp

@@ -21,146 +21,146 @@
 
 
 IGL_INLINE void igl::streamlines_init(
-        const Eigen::MatrixXd V,
-        const Eigen::MatrixXi F,
-        const Eigen::MatrixXd &temp_field,
-        const bool treat_as_symmetric,
-        StreamlineData &data,
-        StreamlineState &state,
-        double percentage
-){
-    using namespace Eigen;
-    using namespace std;
-
-    igl::edge_topology(V, F, data.E, data.F2E, data.E2F);
-    igl::triangle_triangle_adjacency(F, data.TT);
-
-    // prepare vector field
-    // --------------------------
-    int half_degree = temp_field.cols() / 3;
-    int degree = treat_as_symmetric ? half_degree * 2 : half_degree;
-    data.degree = degree;
-
-    Eigen::MatrixXd FN;
-    Eigen::VectorXi order, inv_order_unused;
-    Eigen::RowVectorXd sorted;
-
-    igl::per_face_normals(V, F, FN);
-    data.field.setZero(F.rows(), degree * 3);
-    for (unsigned i = 0; i < F.rows(); ++i){
-        const Eigen::RowVectorXd &n = FN.row(i);
-        Eigen::RowVectorXd temp(1, degree * 3);
-        if (treat_as_symmetric)
-            temp << temp_field.row(i), -temp_field.row(i);
-        else
-            temp = temp_field.row(i);
-        igl::sort_vectors_ccw(temp, n, order, true, sorted, false, inv_order_unused);
-
-        // project vectors to tangent plane
-        for (int j = 0; j < degree; ++j)
-        {
-            Eigen::RowVector3d pd = sorted.segment(j * 3, 3);
-            pd = (pd - (n.dot(pd)) * n).normalized();
-            data.field.block(i, j * 3, 1, 3) = pd;
-        }
+                                      const Eigen::MatrixXd V,
+                                      const Eigen::MatrixXi F,
+                                      const Eigen::MatrixXd &temp_field,
+                                      const bool treat_as_symmetric,
+                                      StreamlineData &data,
+                                      StreamlineState &state,
+                                      double percentage
+                                      ){
+  using namespace Eigen;
+  using namespace std;
+  
+  igl::edge_topology(V, F, data.E, data.F2E, data.E2F);
+  igl::triangle_triangle_adjacency(F, data.TT);
+  
+  // prepare vector field
+  // --------------------------
+  int half_degree = temp_field.cols() / 3;
+  int degree = treat_as_symmetric ? half_degree * 2 : half_degree;
+  data.degree = degree;
+  
+  Eigen::MatrixXd FN;
+  Eigen::VectorXi order;
+  Eigen::RowVectorXd sorted;
+  
+  igl::per_face_normals(V, F, FN);
+  data.field.setZero(F.rows(), degree * 3);
+  for (unsigned i = 0; i < F.rows(); ++i){
+    const Eigen::RowVectorXd &n = FN.row(i);
+    Eigen::RowVectorXd temp(1, degree * 3);
+    if (treat_as_symmetric)
+      temp << temp_field.row(i), -temp_field.row(i);
+    else
+      temp = temp_field.row(i);
+    igl::sort_vectors_ccw(temp, n, order, sorted);
+    
+    // project vectors to tangent plane
+    for (int j = 0; j < degree; ++j)
+    {
+      Eigen::RowVector3d pd = sorted.segment(j * 3, 3);
+      pd = (pd - (n.dot(pd)) * n).normalized();
+      data.field.block(i, j * 3, 1, 3) = pd;
     }
-    Eigen::VectorXd curl;
-    igl::polyvector_field_matchings(data.field, V, F, false, data.match_ab, data.match_ba, curl);
-
-    // create seeds for tracing
-    // --------------------------
-    Eigen::VectorXi samples;
-    int nsamples;
-
-    nsamples = percentage * F.rows();
-    Eigen::VectorXd r;
-    r.setRandom(nsamples, 1);
-    r = (1 + r.array()) / 2.;
-    samples = (r.array() * F.rows()).cast<int>();
-    data.nsample = nsamples;
-
-    Eigen::MatrixXd BC, BC_sample;
-    igl::barycenter(V, F, BC);
-    igl::slice(BC, samples, 1, BC_sample);
-
-    // initialize state for tracing vector field
-
-    state.start_point = BC_sample.replicate(degree,1);
-    state.end_point = state.start_point;
-
-    state.current_face = samples.replicate(1, degree);
-
-    state.current_direction.setZero(nsamples, degree);
-    for (int i = 0; i < nsamples; ++i)
-        for (int j = 0; j < degree; ++j)
-            state.current_direction(i, j) = j;
-
+  }
+  Eigen::VectorXd curl;
+  igl::polyvector_field_matchings(data.field, V, F, false, treat_as_symmetric, data.match_ab, data.match_ba, curl);
+  
+  // create seeds for tracing
+  // --------------------------
+  Eigen::VectorXi samples;
+  int nsamples;
+  
+  nsamples = percentage * F.rows();
+  Eigen::VectorXd r;
+  r.setRandom(nsamples, 1);
+  r = (1 + r.array()) / 2.;
+  samples = (r.array() * F.rows()).cast<int>();
+  data.nsample = nsamples;
+  
+  Eigen::MatrixXd BC, BC_sample;
+  igl::barycenter(V, F, BC);
+  igl::slice(BC, samples, 1, BC_sample);
+  
+  // initialize state for tracing vector field
+  
+  state.start_point = BC_sample.replicate(degree,1);
+  state.end_point = state.start_point;
+  
+  state.current_face = samples.replicate(1, degree);
+  
+  state.current_direction.setZero(nsamples, degree);
+  for (int i = 0; i < nsamples; ++i)
+    for (int j = 0; j < degree; ++j)
+      state.current_direction(i, j) = j;
+  
 }
 
 IGL_INLINE void igl::streamlines_next(
-        const Eigen::MatrixXd V,
-        const Eigen::MatrixXi F,
-        const StreamlineData & data,
-        StreamlineState & state
-){
-    using namespace Eigen;
-    using namespace std;
-
-    int degree = data.degree;
-    int nsample = data.nsample;
-
-    state.start_point = state.end_point;
-
-    for (int i = 0; i < degree; ++i)
+                                      const Eigen::MatrixXd V,
+                                      const Eigen::MatrixXi F,
+                                      const StreamlineData & data,
+                                      StreamlineState & state
+                                      ){
+  using namespace Eigen;
+  using namespace std;
+  
+  int degree = data.degree;
+  int nsample = data.nsample;
+  
+  state.start_point = state.end_point;
+  
+  for (int i = 0; i < degree; ++i)
+  {
+    for (int j = 0; j < nsample; ++j)
     {
-        for (int j = 0; j < nsample; ++j)
+      int f0 = state.current_face(j,i);
+      if (f0 == -1) // reach boundary
+        continue;
+      int m0 = state.current_direction(j, i);
+      
+      // the starting point of the vector
+      const Eigen::RowVector3d &p = state.start_point.row(j + nsample * i);
+      // the direction where we are trying to go
+      const Eigen::RowVector3d &r = data.field.block(f0, 3 * m0, 1, 3);
+      
+      
+      // new state,
+      int f1, m1;
+      
+      for (int k = 0; k < 3; ++k)
+      {
+        f1 = data.TT(f0, k);
+        
+        // edge vertices
+        const Eigen::RowVector3d &q = V.row(F(f0, k));
+        const Eigen::RowVector3d &qs = V.row(F(f0, (k + 1) % 3));
+        // edge direction
+        Eigen::RowVector3d s = qs - q;
+        
+        double u;
+        double t;
+        if (igl::segments_intersect(p, r, q, s, t, u))
         {
-            int f0 = state.current_face(j,i);
-            if (f0 == -1) // reach boundary
-                continue;
-            int m0 = state.current_direction(j, i);
-
-            // the starting point of the vector
-            const Eigen::RowVector3d &p = state.start_point.row(j + nsample * i);
-            // the direction where we are trying to go
-            const Eigen::RowVector3d &r = data.field.block(f0, 3 * m0, 1, 3);
-
-
-            // new state,
-            int f1, m1;
-
-            for (int k = 0; k < 3; ++k)
-            {
-                f1 = data.TT(f0, k);
-
-                // edge vertices
-                const Eigen::RowVector3d &q = V.row(F(f0, k));
-                const Eigen::RowVector3d &qs = V.row(F(f0, (k + 1) % 3));
-                // edge direction
-                Eigen::RowVector3d s = qs - q;
-
-                double u;
-                double t;
-                if (igl::segments_intersect(p, r, q, s, t, u))
-                {
-                    // point on next face
-                    state.end_point.row(j + nsample * i) = p + t * r;
-                    state.current_face(j,i) = f1;
-
-                    // matching direction on next face
-                    int e1 = data.F2E(f0, k);
-                    if (data.E2F(e1, 0) == f0)
-                        m1 = data.match_ab(e1, m0);
-                    else
-                        m1 = data.match_ba(e1, m0);
-
-                    state.current_direction(j, i) = m1;
-                    break;
-                }
-
-            }
-
-
+          // point on next face
+          state.end_point.row(j + nsample * i) = p + t * r;
+          state.current_face(j,i) = f1;
+          
+          // matching direction on next face
+          int e1 = data.F2E(f0, k);
+          if (data.E2F(e1, 0) == f0)
+            m1 = data.match_ab(e1, m0);
+          else
+            m1 = data.match_ba(e1, m0);
+          
+          state.current_direction(j, i) = m1;
+          break;
         }
+        
+      }
+      
+      
     }
+  }
 }

include/igl/triangle_triangle_adjacency.cpp

@@ -220,4 +220,5 @@ template void igl::triangle_triangle_adjacency<Eigen::Matrix<int, -1, -1, 0, -1,
 template void igl::triangle_triangle_adjacency<Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, long, long, long>(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> > const&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > > const&, bool, std::vector<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >, std::allocator<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > > > >&, std::vector<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >, std::allocator<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > > > >&);
 template void igl::triangle_triangle_adjacency<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, unsigned long, int, int>(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, std::vector<std::vector<unsigned long, std::allocator<unsigned long> >, std::allocator<std::vector<unsigned long, std::allocator<unsigned long> > > > const&, bool, std::vector<std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >, std::allocator<std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > > >&, std::vector<std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >, std::allocator<std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > > >&);
 template void igl::triangle_triangle_adjacency<Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
+template void igl::triangle_triangle_adjacency<Eigen::Matrix<int, -1, -1, 0, -1, -1>, long, long>(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >, std::allocator<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > > > >&, std::vector<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >, std::allocator<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > > > >&);
 #endif

include/igl/unique_edge_map.cpp

@@ -55,6 +55,7 @@ template void igl::unique_edge_map<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen:
 
 #ifdef WIN32
 template void __cdecl igl::unique_edge_map<class Eigen::Matrix<int, -1, 3, 0, -1, 3>, class Eigen::Matrix<int, -1, 2, 0, -1, 2>, class Eigen::Matrix<int, -1, 2, 0, -1, 2>, class Eigen::Matrix<__int64, -1, 1, 0, -1, 1>, __int64>(class Eigen::PlainObjectBase<class Eigen::Matrix<int, -1, 3, 0, -1, 3> > const &, class Eigen::PlainObjectBase<class Eigen::Matrix<int, -1, 2, 0, -1, 2> > &, class Eigen::PlainObjectBase<class Eigen::Matrix<int, -1, 2, 0, -1, 2> > &, class Eigen::PlainObjectBase<class Eigen::Matrix<__int64, -1, 1, 0, -1, 1> > &, class std::vector<class std::vector<__int64, class std::allocator<__int64> >, class std::allocator<class std::vector<__int64, class std::allocator<__int64> > > > &);
+template void __cdecl igl::unique_edge_map<class Eigen::Matrix<int,-1,-1,0,-1,-1>,class Eigen::Matrix<int,-1,2,0,-1,2>,class Eigen::Matrix<int,-1,2,0,-1,2>,class Eigen::Matrix<__int64,-1,1,0,-1,1>,__int64>(class Eigen::PlainObjectBase<class Eigen::Matrix<int,-1,-1,0,-1,-1> > const &,class Eigen::PlainObjectBase<class Eigen::Matrix<int,-1,2,0,-1,2> > &,class Eigen::PlainObjectBase<class Eigen::Matrix<int,-1,2,0,-1,2> > &,class Eigen::PlainObjectBase<class Eigen::Matrix<__int64,-1,1,0,-1,1> > &,class std::vector<class std::vector<__int64,class std::allocator<__int64> >,class std::allocator<class std::vector<__int64,class std::allocator<__int64> > > > &);
 #endif
 
 #endif 

optional/CMakeLists.txt

@@ -1,6 +1,10 @@
 cmake_minimum_required(VERSION 2.8.12)
 project(libigl)
 
+set (CMAKE_MODULE_PATH
+#${CMAKE_MODULE_PATH}
+"${PROJECT_SOURCE_DIR}/../shared/cmake")
+
 ### Compilation flags: adapt to your needs ###
 if(MSVC)
   set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /MP /bigobj /w") ### Enable parallel compilation
@@ -31,7 +35,7 @@ option(LIBIGL_WITH_MATLAB           "Use Matlab"         "${MATLAB_FOUND}")
 find_package(MOSEK  QUIET)
 option(LIBIGL_WITH_MOSEK            "Use MOSEK"          "${MOSEK_FOUND}")
 ### Nanogui is off by default because it has many dependencies and generates
-### many issues 
+### many issues
 option(LIBIGL_WITH_NANOGUI          "Use Nanogui menu"   OFF)
 option(LIBIGL_WITH_OPENGL           "Use OpenGL"         ON)
 option(LIBIGL_WITH_PNG              "Use PNG"            ON)

shared/cmake/FindMATLAB.cmake

@@ -137,7 +137,7 @@ ELSE(WIN32)
     IF((NOT DEFINED MATLAB_ROOT) OR ("${MATLAB_ROOT}" STREQUAL ""))
 
     # Search for a version of Matlab available, starting from the most modern one to older versions
-      FOREACH(MATVER "R2015b" "R2015a" "R2014b" "R2014a" "R2014a" "R2013b" "R2013a" "R2012b" "R2012a" "R2011b" "R2011a" "R2010b" "R2010a" "R2009b" "R2009a" "R2008b")
+      FOREACH(MATVER "R2016a" "R2015b" "R2015a" "R2014b" "R2014a" "R2014a" "R2013b" "R2013a" "R2012b" "R2012a" "R2011b" "R2011a" "R2010b" "R2010a" "R2009b" "R2009a" "R2008b")
         IF((NOT DEFINED MATLAB_ROOT) OR ("${MATLAB_ROOT}" STREQUAL ""))
           IF(EXISTS /Applications/MATLAB_${MATVER}.app)
             SET(MATLAB_ROOT /Applications/MATLAB_${MATVER}.app)

tutorial/508_ConjugateField/main.cpp

@@ -33,8 +33,6 @@ Eigen::MatrixXd conjugate_pvf;
 Eigen::VectorXd conjugacy_s;
 Eigen::VectorXd conjugacy_c;
 
-igl::ConjugateFFSolverData<Eigen::MatrixXd, Eigen::MatrixXi> *csdata;
-
 
 bool key_down(igl::viewer::Viewer& viewer, unsigned char key, int modifier)
 {
@@ -143,7 +141,7 @@ int main(int argc, char *argv[])
   igl::n_polyvector(V, F, b, bc, smooth_pvf);
 
   // Initialize conjugate field with smooth field
-  csdata = new igl::ConjugateFFSolverData<Eigen::MatrixXd,Eigen::MatrixXi>(V,F);
+  igl::ConjugateFFSolverData<Eigen::MatrixXd, Eigen::MatrixXi> csdata(V,F);
   conjugate_pvf = smooth_pvf;
 
 
@@ -153,12 +151,11 @@ int main(int argc, char *argv[])
   double lambdaInit = 100;
   double lambdaMultFactor = 1.01;
   bool doHardConstraints = true;
-  double lambdaOut;
   VectorXi isConstrained = VectorXi::Constant(F.rows(),0);
   for (unsigned i=0; i<b.size(); ++i)
     isConstrained(b(i)) = 1;
 
-  igl::conjugate_frame_fields(*csdata, isConstrained, conjugate_pvf, conjugate_pvf, conjIter, lambdaOrtho, lambdaInit, lambdaMultFactor, doHardConstraints, &lambdaOut);
+  double lambdaOut = igl::conjugate_frame_fields(csdata, isConstrained, conjugate_pvf, conjugate_pvf, conjIter, lambdaOrtho, lambdaInit, lambdaMultFactor, doHardConstraints);
 
   // local representations of field vectors
   Eigen::Matrix<double, Eigen::Dynamic, 2> pvU, pvV;
@@ -168,13 +165,13 @@ int main(int argc, char *argv[])
   const Eigen::MatrixXd &Vs = smooth_pvf.rightCols(3);
   pvU << igl::dot_row(Us,B1), igl::dot_row(Us,B2);
   pvV << igl::dot_row(Vs,B1), igl::dot_row(Vs,B2);
-  csdata->evaluateConjugacy(pvU, pvV, conjugacy_s);
+  csdata.evaluateConjugacy(pvU, pvV, conjugacy_s);
   //conjugate
   const Eigen::MatrixXd &Uc = conjugate_pvf.leftCols(3);
   const Eigen::MatrixXd &Vc = conjugate_pvf.rightCols(3);
   pvU << igl::dot_row(Uc,B1), igl::dot_row(Uc,B2);
   pvV << igl::dot_row(Vc,B1), igl::dot_row(Vc,B2);
-  csdata->evaluateConjugacy(pvU, pvV, conjugacy_c);
+  csdata.evaluateConjugacy(pvU, pvV, conjugacy_c);
   // Launch the viewer
   igl::viewer::Viewer viewer;
   viewer.core.invert_normals = true;

tutorial/510_Integrable/main.cpp

@@ -595,7 +595,7 @@ bool key_down(igl::viewer::Viewer& viewer, unsigned char key, int modifier)
     // Compute curl_minimizing matchings and curl
     printf("--Matchings and curl--\n");
     Eigen::MatrixXi match_ab, match_ba;  // matchings across interior edges
-    double avgCurl = igl::polyvector_field_matchings(two_pv, V, F, true, match_ab, match_ba, curl);
+    double avgCurl = igl::polyvector_field_matchings(two_pv, V, F, true, true, match_ab, match_ba, curl);
     double maxCurl = curl.maxCoeff();
     printf("curl -- max: %.5g, avg: %.5g\n", maxCurl,  avgCurl);
     // Compute singularities
@@ -662,7 +662,7 @@ int main(int argc, char *argv[])
   // Compute curl_minimizing matchings and curl
   Eigen::MatrixXi match_ab, match_ba;  // matchings across interior edges
   printf("--Matchings and curl--\n");
-  double avgCurl = igl::polyvector_field_matchings(two_pv_ori, V, F, true, match_ab, match_ba, curl_ori);
+  double avgCurl = igl::polyvector_field_matchings(two_pv_ori, V, F, true, true, match_ab, match_ba, curl_ori);
   double maxCurl = curl_ori.maxCoeff();
   printf("original curl -- max: %.5g, avg: %.5g\n", maxCurl,  avgCurl);
 

tutorial/511_PolyVectorFieldGeneral/CMakeLists.txt

@@ -0,0 +1,8 @@
+cmake_minimum_required(VERSION 2.8.12)
+project(511_PolyVectorFieldGeneral)
+
+add_executable(${PROJECT_NAME}_bin
+  main.cpp)
+target_include_directories(${PROJECT_NAME}_bin PRIVATE ${LIBIGL_INCLUDE_DIRS})
+target_compile_definitions(${PROJECT_NAME}_bin PRIVATE ${LIBIGL_DEFINITIONS})
+target_link_libraries(${PROJECT_NAME}_bin ${LIBIGL_LIBRARIES} ${LIBIGL_EXTRA_LIBRARIES})

tutorial/511_PolyVectorFieldGeneral/main.cpp

@@ -0,0 +1,183 @@
+#include <igl/readOBJ.h>
+#include <igl/readDMAT.h>
+#include <igl/viewer/Viewer.h>
+#include <igl/barycenter.h>
+#include <igl/avg_edge_length.h>
+#include <vector>
+#include <igl/n_polyvector_general.h>
+#include <igl/n_polyvector.h>
+#include <igl/local_basis.h>
+#include <igl/writeOFF.h>
+#include <stdlib.h>
+#include <igl/jet.h>
+#include <fstream>
+#include <iostream>
+// Input mesh
+Eigen::MatrixXd V;
+Eigen::MatrixXi F;
+
+// Per face bases
+Eigen::MatrixXd B1,B2,B3;
+
+// Face barycenters
+Eigen::MatrixXd B;
+
+// Scale for visualizing the fields
+double global_scale;
+
+// Random length factor
+double rand_factor = 5;
+
+Eigen::VectorXi samples;
+
+void readSamples(const std::string &fname, Eigen::VectorXi &samples)
+{
+  int numSamples;
+  FILE *fp = fopen(fname.c_str(),"r");
+  if (fscanf(fp, "%d", &numSamples)!=1)
+  {
+    fclose(fp);
+    return;
+  }
+  samples.resize(numSamples,1);
+  int vali;
+  for (int i =0; i<numSamples; ++i)
+  {
+    if (fscanf(fp, "%d", &vali)!=1 || vali<0)
+    {
+      fclose(fp);
+      samples.resize(0,1);
+      return;
+    }
+    samples[i]=vali;
+  }
+  fclose(fp);
+  
+}
+
+// Create a random set of tangent vectors
+Eigen::VectorXd random_constraints(const
+                                   Eigen::VectorXd& b1, const
+                                   Eigen::VectorXd& b2, int n)
+{
+  Eigen::VectorXd r(n*3);
+  for (unsigned i=0; i<n;++i)
+  {
+    double a = (double(rand())/RAND_MAX)*2*M_PI;
+    double s = 1 + ((double(rand())/RAND_MAX)) * rand_factor;
+    Eigen::Vector3d t = s * (cos(a) * b1 + sin(a) * b2);
+    r.block(i*3,0,3,1) = t;
+  }
+  return r;
+}
+
+bool key_down(igl::viewer::Viewer& viewer, unsigned char key, int modifier)
+{
+  using namespace std;
+  using namespace Eigen;
+
+  if (key <'1' || key >'9')
+    return false;
+
+  viewer.data.lines.resize(0,9);
+
+  int num = key  - '0';
+
+  // Interpolate
+  std::cerr << "Interpolating " << num  << "-PolyVector field" << std::endl;
+
+  VectorXi b(3);
+  b << 1511, 603, 506;
+
+  int numConstraintsToGenerate;
+  // if it's not a 2-PV or a 1-PV, include a line direction (2 opposite vectors)
+  // in the field
+  if (num>=5)
+    numConstraintsToGenerate  = num-2;
+  else
+    if (num>=3)
+      numConstraintsToGenerate  = num-1;
+    else
+      numConstraintsToGenerate  = num;
+
+
+  MatrixXd bc(b.size(),numConstraintsToGenerate*3);
+  for (unsigned i=0; i<b.size(); ++i)
+  {
+    VectorXd t = random_constraints(B1.row(b(i)),B2.row(b(i)),numConstraintsToGenerate);
+    bc.row(i) = t;
+  }
+  VectorXi rootsIndex(num);
+  for (int i =0; i<numConstraintsToGenerate; ++i)
+    rootsIndex[i] = i+1;
+  if (num>=5)
+    rootsIndex[num-2] = -2;
+    if (num>=3)
+      rootsIndex[num-1] = -1;
+
+  // Interpolated PolyVector field
+  Eigen::MatrixXd pvf;
+  igl::n_polyvector_general(V, F, b, bc, rootsIndex, pvf);
+
+  ofstream ofs;
+  ofs.open("pvf.txt", ofstream::out);
+  ofs<<pvf;
+  ofs.close();
+  igl::writeOFF("pvf.off",V,F);
+  
+  // Highlight in red the constrained faces
+  MatrixXd C = MatrixXd::Constant(F.rows(),3,1);
+  for (unsigned i=0; i<b.size();++i)
+    C.row(b(i)) << 1, 0, 0;
+  viewer.data.set_colors(C);
+
+  for (int n=0; n<num; ++n)
+  {
+//    const MatrixXd &VF = pvf.block(0,n*3,F.rows(),3);
+    MatrixXd VF = MatrixXd::Zero(F.rows(),3);
+    for (unsigned i=0; i<b.size(); ++i)
+      VF.row(b[i]) = pvf.block(b[i],n*3,1,3);
+    
+    for (int i=0; i<samples.rows(); ++i)
+      VF.row(samples[i]) = pvf.block(samples[i],n*3,1,3);
+    
+    VectorXd c = VF.rowwise().norm();
+    MatrixXd C2;
+    igl::jet(c,1,1+rand_factor,C2);
+    // viewer.data.add_edges(B - global_scale*VF, B + global_scale*VF , C2);
+    viewer.data.add_edges(B, B + global_scale*VF , C2);
+  }
+
+
+  return false;
+}
+
+int main(int argc, char *argv[])
+{
+  using namespace Eigen;
+  using namespace std;
+  // Load a mesh in OBJ format
+  igl::readOBJ(TUTORIAL_SHARED_PATH "/snail.obj", V, F);
+  readSamples(TUTORIAL_SHARED_PATH "/snail.samples.0.2", samples);
+
+  // Compute local basis for faces
+  igl::local_basis(V,F,B1,B2,B3);
+
+  // Compute face barycenters
+  igl::barycenter(V, F, B);
+
+  // Compute scale for visualizing fields
+  global_scale =  .1*igl::avg_edge_length(V, F);
+
+  // Make the example deterministic
+  srand(0);
+
+  igl::viewer::Viewer viewer;
+  viewer.data.set_mesh(V, F);
+  viewer.callback_key_down = &key_down;
+  viewer.core.show_lines = false;
+  key_down(viewer,'3',0);
+  std::cerr << " *** Press keys 1-9 to select number of vectors per point. ***" << std::endl;
+  
+  viewer.launch();
+}

tutorial/CMakeLists.txt

@@ -137,6 +137,7 @@ if(TUTORIALS_CHAPTER5)
   add_subdirectory("508_ConjugateField")
   add_subdirectory("509_Planarization")
   add_subdirectory("510_Integrable")
+  add_subdirectory("511_PolyVectorFieldGeneral")
 endif()
 
 # Chapter 6

tutorial/images/511_PolyVectorFieldGeneral.png.REMOVED.git-id

@@ -0,0 +1 @@
+9a6cb4e375e2b9a94155e147f94dfe3684c23265

tutorial/tutorial.md.REMOVED.git-id

@@ -1 +1 @@
-d0e101aea39668e374ff776b29cdd79667b7d8bd
+0e9dc3a10233852a045d53b248481da495b37e4c