Merge branch 'master' of https://github.com/olkido/libigl

Conflicts:
	tutorial/cmake/FindLIBCOMISO.cmake

Former-commit-id: 199b4aa66be125041238c0391b96e7f3450d958b

include/igl/conjugate_frame_fields.cpp

@@ -7,26 +7,79 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 
 #include <igl/conjugate_frame_fields.h>
-#include <igl/edgetopology.h>
+#include <igl/edge_topology.h>
 #include <igl/local_basis.h>
 #include <igl/nchoosek.h>
 #include <igl/sparse.h>
 #include <igl/speye.h>
 #include <igl/slice.h>
 #include <igl/polyroots.h>
-#include <igl/add_barycenter.h>
+#include <igl/colon.h>
+#include <igl/false_barycentric_subdivision.h>
 #include <igl/principal_curvature.h>
 #include <Eigen/Sparse>
 
 #include <iostream>
 
 namespace igl {
+
+  template <typename DerivedV, typename DerivedF>
+  class ConjugateFFSolverData
+  {
+    public:
+      const Eigen::PlainObjectBase<DerivedV> &V; int numV;
+      const Eigen::PlainObjectBase<DerivedF> &F; int numF;
+
+      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;
+
+      Eigen::PlainObjectBase<DerivedV> B1, B2, FN;
+
+
+      Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic,1> kmin, kmax;
+      Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic,2> dmin, dmax;
+      Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic,3> dmin3, dmax3;
+
+      Eigen::VectorXd nonPlanarityMeasure;
+      Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > planarityWeight;
+
+      //conjugacy matrix
+      std::vector<Eigen::Matrix<typename DerivedV::Scalar, 4,4> > H;
+
+      //conjugacy matrix eigenvectors and (scaled) eigenvalues
+      std::vector<Eigen::Matrix<typename DerivedV::Scalar, 4,4> > UH;
+      std::vector<Eigen::Matrix<typename DerivedV::Scalar, 4,1> > s;
+
+      //laplacians
+      Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar>> DDA, DDB;
+
+  private:
+    IGL_INLINE void computeCurvatureAndPrincipals();
+    IGL_INLINE void precomputeConjugacyStuff();
+    IGL_INLINE void computeLaplacians();
+    IGL_INLINE void computek();
+    IGL_INLINE void computeCoefficientLaplacian(int n, Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &D);
+    
+    IGL_INLINE void precomputeInteriorEdges();
+
+public:
+      IGL_INLINE ConjugateFFSolverData(const Eigen::PlainObjectBase<DerivedV> &_V,
+                                   const Eigen::PlainObjectBase<DerivedF> &_F);
+  };
+
   template <typename DerivedV, typename DerivedF, typename DerivedO>
   class ConjugateFFSolver
   {
   public:
-    IGL_INLINE ConjugateFFSolver(const Eigen::PlainObjectBase<DerivedV> &_V,
-                                 const Eigen::PlainObjectBase<DerivedF> &_F,
+    IGL_INLINE ConjugateFFSolver(const ConjugateFFSolverData<DerivedV, DerivedF> &_data,
                                  int _maxIter = 50,
                                  const typename DerivedV::Scalar &_lambdaOrtho = .1,
                                  const typename DerivedV::Scalar &_lambdaInit = 100,
@@ -34,72 +87,35 @@ namespace igl {
                                  bool _doHardConstraints = true);
     IGL_INLINE bool solve(const Eigen::VectorXi &isConstrained,
                           const Eigen::PlainObjectBase<DerivedO> &initialSolution,
-                          Eigen::PlainObjectBase<DerivedO> &output);
-    
+                          Eigen::PlainObjectBase<DerivedO> &output,
+                          typename DerivedV::Scalar *lambdaOut = NULL);
+
   private:
-    const Eigen::PlainObjectBase<DerivedV> &V; int numV;
-    const Eigen::PlainObjectBase<DerivedF> &F; int numF;
-    
-    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;
-    
-    Eigen::PlainObjectBase<DerivedV> B1, B2, FN;
-    
-    
-    Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic,1> kmin, kmax;
-    Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic,2> dmin, dmax;
-    Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic,3> dmin3, dmax3;
-    
-    Eigen::VectorXd nonPlanarityMeasure;
-    Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > planarityWeight;
-    
-    //conjugacy matrix
-    std::vector<Eigen::Matrix<typename DerivedV::Scalar, 4,4> > H;
-    
-    //conjugacy matrix eigenvectors and (scaled) eigenvalues
-    std::vector<Eigen::Matrix<typename DerivedV::Scalar, 4,4> > UH;
-    std::vector<Eigen::Matrix<typename DerivedV::Scalar, 4,1> > s;
-    
-    //laplacians
-    Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar>> DDA, DDB;
-    
+
+    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 computeCurvatureAndPrincipals();
+
+
+
     IGL_INLINE void evaluateConjugacy(Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 1> &conjValues);
-    
-    IGL_INLINE void precomputeConjugacyStuff();
-    IGL_INLINE void computeLaplacians();
-    IGL_INLINE void computek();
-    IGL_INLINE void computeCoefficientLaplacian(int n, Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &D);
-    
-    IGL_INLINE void precomputeInteriorEdges();
-    
-    
+
+
+
     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);
@@ -114,70 +130,57 @@ namespace igl {
   };
 }
 
-template <typename DerivedV, typename DerivedF, typename DerivedO>
-IGL_INLINE igl::ConjugateFFSolver<DerivedV, DerivedF, DerivedO>::
-ConjugateFFSolver(const Eigen::PlainObjectBase<DerivedV> &_V,
-                  const Eigen::PlainObjectBase<DerivedF> &_F,
-                  int _maxIter,
-                  const typename DerivedV::Scalar &_lambdaOrtho,
-                  const typename DerivedV::Scalar &_lambdaInit,
-                  const typename DerivedV::Scalar &_lambdaMultFactor,
-                  bool _doHardConstraints):
+//Implementation
+/***************************** Data ***********************************/
+
+template <typename DerivedV, typename DerivedF>
+IGL_INLINE igl::ConjugateFFSolverData<DerivedV, DerivedF>::
+ConjugateFFSolverData(const Eigen::PlainObjectBase<DerivedV> &_V,
+                  const Eigen::PlainObjectBase<DerivedF> &_F):
 V(_V),
 numV(_V.rows()),
 F(_F),
-numF(_F.rows()),
-lambdaOrtho(_lambdaOrtho),
-lambdaInit(_lambdaInit),
-maxIter(_maxIter),
-lambdaMultFactor(_lambdaMultFactor),
-doHardConstraints(_doHardConstraints)
+numF(_F.rows())
 {
-  
-  igl::edgetopology(V,F,EV,F2E,E2F);
+  igl::edge_topology(V,F,EV,F2E,E2F);
   numE = EV.rows();
-  
+
   precomputeInteriorEdges();
-  
+
   igl::local_basis(V,F,B1,B2,FN);
-  
+
   computek();
-  
+
   computeLaplacians();
-  
+
   computeCurvatureAndPrincipals();
   precomputeConjugacyStuff();
-  
-  Acoeff.resize(numF,1);
-  Bcoeff.resize(numF,1);
-  pvU.setZero(numF, 2);
-  pvV.setZero(numF, 2);
 
 };
 
 
-template <typename DerivedV, typename DerivedF, typename DerivedO>
-IGL_INLINE void igl::ConjugateFFSolver<DerivedV, DerivedF, DerivedO>::computeCurvatureAndPrincipals()
+template <typename DerivedV, typename DerivedF>
+IGL_INLINE void igl::ConjugateFFSolverData<DerivedV, DerivedF>::computeCurvatureAndPrincipals()
 {
   Eigen::MatrixXd VCBary;
   Eigen::MatrixXi FCBary;
-  
+
   VCBary.setZero(numV+numF,3);
   FCBary.setZero(3*numF,3);
-  igl::add_barycenter(V, F, VCBary, FCBary);
-  
+  igl::false_barycentric_subdivision(V, F, VCBary, FCBary);
+
   Eigen::MatrixXd dmax3_,dmin3_;
   igl::principal_curvature(VCBary, FCBary, dmax3_, dmin3_, kmax, kmin, 5,true);
-  
+
   dmax3 = dmax3_.bottomRows(numF);
   dmin3 = dmin3_.bottomRows(numF);
-  
+
   kmax = kmax.bottomRows(numF);
   kmin = kmin.bottomRows(numF);
-  
+
   //  kmax = dmax3.rowwise().norm();
   //  kmin = dmin3.rowwise().norm();
-  
+
   dmin3.rowwise().normalize();
   dmax3.rowwise().normalize();
   dmax.setZero(numF,2);
@@ -203,30 +206,30 @@ IGL_INLINE void igl::ConjugateFFSolver<DerivedV, DerivedF, DerivedO>::computeCur
     dmax.row(i).normalize();
     dmin.row(i) << dmin3.row(i).dot(B1.row(i)), dmin3.row(i).dot(B2.row(i));
     dmin.row(i).normalize();
-    
+
   }
-  
+
   nonPlanarityMeasure = kmax.cwiseAbs().array()*kmin.cwiseAbs().array();
   typename DerivedV::Scalar minP = nonPlanarityMeasure.minCoeff();
   typename DerivedV::Scalar maxP = nonPlanarityMeasure.maxCoeff();
   nonPlanarityMeasure = (nonPlanarityMeasure.array()-minP)/(maxP-minP);
   Eigen::VectorXi I = igl::colon<typename DerivedF::Scalar>(0, numF-1);
   igl::sparse(I, I, nonPlanarityMeasure, numF, numF, planarityWeight);
-  
+
 }
 
-template <typename DerivedV, typename DerivedF, typename DerivedO>
-IGL_INLINE void igl::ConjugateFFSolver<DerivedV, DerivedF, DerivedO>::precomputeConjugacyStuff()
+template <typename DerivedV, typename DerivedF>
+IGL_INLINE void igl::ConjugateFFSolverData<DerivedV, DerivedF>::precomputeConjugacyStuff()
 {
   H.resize(numF);
   UH.resize(numF);
   s.resize(numF);
-  
+
   for (int i = 0; i<numF; ++i)
   {
     //compute conjugacy matrix
     typename DerivedV::Scalar e1x = dmin(i,0), e1y = dmin(i,1), e2x = dmax(i,0), e2y = dmax(i,1), k1 = kmin[i], k2 = kmax[i];
-    
+
     H[i]<<
     0,          0, k1*e1x*e1x, k1*e1x*e1y,
     0,          0, k1*e1x*e1y, k1*e1y*e1y,
@@ -234,35 +237,35 @@ IGL_INLINE void igl::ConjugateFFSolver<DerivedV, DerivedF, DerivedO>::precompute
     k2*e2x*e2y, k2*e2y*e2y,          0,          0;
     Eigen::Matrix<typename DerivedV::Scalar, 4, 4> Ht = H[i].transpose();
     H[i] = .5*(H[i]+Ht);
-    
+
     Eigen::EigenSolver<Eigen::Matrix<typename DerivedV::Scalar, 4, 4> > es(H[i]);
     s[i] = es.eigenvalues().real();//ok to do this because H symmetric
     //scale
     s[i] = s[i]/(s[i].cwiseAbs().minCoeff());
     UH[i] = es.eigenvectors().real();
-    
-    
+
+
   }
 }
 
 
-template <typename DerivedV, typename DerivedF, typename DerivedO>
-IGL_INLINE void igl::ConjugateFFSolver<DerivedV, DerivedF, DerivedO>::computeLaplacians()
+template <typename DerivedV, typename DerivedF>
+IGL_INLINE void igl::ConjugateFFSolverData<DerivedV, DerivedF>::computeLaplacians()
 {
   computeCoefficientLaplacian(2, DDA);
-  
+
   computeCoefficientLaplacian(4, DDB);
 }
 
-template<typename DerivedV, typename DerivedF, typename DerivedO>
-IGL_INLINE void igl::ConjugateFFSolver<DerivedV, DerivedF, DerivedO>::
+template<typename DerivedV, typename DerivedF>
+IGL_INLINE void igl::ConjugateFFSolverData<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)))
@@ -273,7 +276,7 @@ precomputeInteriorEdges()
       numInteriorEdges++;
     }
   }
-  
+
   E2F_int.resize(numInteriorEdges, 2);
   indInteriorToFull.setZero(numInteriorEdges,1);
   int ii = 0;
@@ -285,18 +288,165 @@ precomputeInteriorEdges()
     indInteriorToFull[ii] = k;
     ii++;
   }
-  
+
+}
+
+
+
+template<typename DerivedV, typename DerivedF>
+IGL_INLINE void igl::ConjugateFFSolverData<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());
+
+
 }
 
+template<typename DerivedV, typename DerivedF>
+IGL_INLINE void igl::ConjugateFFSolverData<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();
+
+      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();
+
+      // 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();
+
+      // 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();
+
+      K[eid] = ktemp;
+    }
+  }
+
+}
+
+
+/***************************** Solver ***********************************/
+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,
+                  bool _doHardConstraints):
+data(_data),
+lambdaOrtho(_lambdaOrtho),
+lambdaInit(_lambdaInit),
+maxIter(_maxIter),
+lambdaMultFactor(_lambdaMultFactor),
+doHardConstraints(_doHardConstraints)
+{
+  Acoeff.resize(data.numF,1);
+  Bcoeff.resize(data.numF,1);
+  pvU.setZero(data.numF, 2);
+  pvV.setZero(data.numF, 2);
+};
+
+
+
 template<typename DerivedV, typename DerivedF, typename DerivedO>
 IGL_INLINE void igl::ConjugateFFSolver<DerivedV, DerivedF, DerivedO>::
 evaluateConjugacy(Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 1> &conjValues)
 {
-  conjValues.resize(numF,1);
-  for (int j =0; j<numF; ++j)
+  conjValues.resize(data.numF,1);
+  for (int j =0; j<data.numF; ++j)
   {
     Eigen::Matrix<typename DerivedV::Scalar, 4, 1> x; x<<pvU.row(j).transpose(), pvV.row(j).transpose();
-    conjValues[j] = x.transpose()*H[j]*x;
+    conjValues[j] = x.transpose()*data.H[j]*x;
   }
 }
 
@@ -315,7 +465,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 ;
@@ -324,7 +474,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;
-  
+
 }
 
 
@@ -332,34 +482,34 @@ template<typename DerivedV, typename DerivedF, typename DerivedO>
 IGL_INLINE void igl::ConjugateFFSolver<DerivedV, DerivedF, DerivedO>::
 localStep()
 {
-  for (int j =0; j<numF; ++j)
+  for (int j =0; j<data.numF; ++j)
   {
     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 = UH[j].transpose()*xproj;
+    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(s[j], z, 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)
     {
       if (fabs(imag(roots[i]))>1e-10)
         continue;
-      Eigen::Matrix<typename DerivedV::Scalar, 4, 4> D = ((Eigen::Matrix<typename DerivedV::Scalar, 4, 1>::Ones()+real(roots(i))*s[j]).array().inverse()).matrix().asDiagonal();
-      Eigen::Matrix<typename DerivedV::Scalar, 4, 1> candidate = UH[j]*D*z;
+      Eigen::Matrix<typename DerivedV::Scalar, 4, 4> D = ((Eigen::Matrix<typename DerivedV::Scalar, 4, 1>::Ones()+real(roots(i))*data.s[j]).array().inverse()).matrix().asDiagonal();
+      Eigen::Matrix<typename DerivedV::Scalar, 4, 1> candidate = data.UH[j]*D*z;
       typename DerivedV::Scalar dist = (candidate-xproj).norm();
       if (dist<minDist)
       {
         minDist = dist;
         x = candidate;
       }
-      
+
     }
-    
+
     pvU.row(j) << x(0),x(1);
     pvV.row(j) << x(2),x(3);
   }
@@ -370,7 +520,7 @@ template<typename DerivedV, typename DerivedF, typename DerivedO>
 IGL_INLINE void igl::ConjugateFFSolver<DerivedV, DerivedF, DerivedO>::
 setCoefficientsFromField()
 {
-  for (int i = 0; i <numF; ++i)
+  for (int i = 0; i <data.numF; ++i)
   {
     std::complex<typename DerivedV::Scalar> u(pvU(i,0),pvU(i,1));
     std::complex<typename DerivedV::Scalar> v(pvV(i,0),pvV(i,1));
@@ -387,15 +537,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(numF, numF, I);
-  Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > QA = DDA+lambda*planarityWeight+lambdaOrtho*I;
-  Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > fA = (-2*lambda*planarityWeight*Acoeff).sparseView();
-  
-  Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > QB = DDB+lambda*planarityWeight;
-  Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > fB = (-2*lambda*planarityWeight*Bcoeff).sparseView();
-  
+  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);
@@ -403,13 +553,13 @@ globalStep(const Eigen::Matrix<int, Eigen::Dynamic, 1>  &isConstrained,
   }
   else
   {
-    Eigen::Matrix<int, Eigen::Dynamic, 1>isknown_; isknown_.setZero(numF,1);
+    Eigen::Matrix<int, Eigen::Dynamic, 1>isknown_; isknown_.setZero(data.numF,1);
     Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> xknown_; xknown_.setZero(0,1);
     minQuadWithKnownMini(QA, fA, isknown_, xknown_, Acoeff);
     minQuadWithKnownMini(QB, fB, isknown_, xknown_, Bcoeff);
   }
   setFieldFromCoefficients();
-  
+
 }
 
 
@@ -417,16 +567,16 @@ template<typename DerivedV, typename DerivedF, typename DerivedO>
 IGL_INLINE void igl::ConjugateFFSolver<DerivedV, DerivedF, DerivedO>::
 setFieldFromCoefficients()
 {
-  for (int i = 0; i <numF; ++i)
+  for (int i = 0; i <data.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(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;
@@ -443,7 +593,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>
@@ -455,11 +605,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])
@@ -472,21 +622,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)
@@ -500,7 +650,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)
@@ -508,7 +658,7 @@ minQuadWithKnownMini(const Eigen::SparseMatrix<std::complex<typename DerivedV::S
       x[i] = xknown[indk++];
     else
       x[i] = b.coeff(indu++,0);
-  
+
 }
 
 
@@ -516,18 +666,19 @@ template<typename DerivedV, typename DerivedF, typename DerivedO>
 IGL_INLINE bool igl::ConjugateFFSolver<DerivedV, DerivedF, DerivedO>::
 solve(const Eigen::VectorXi &isConstrained,
       const Eigen::PlainObjectBase<DerivedO> &initialSolution,
-      Eigen::PlainObjectBase<DerivedO> &output)
+      Eigen::PlainObjectBase<DerivedO> &output,
+      typename DerivedV::Scalar *lambdaOut)
 {
   int numConstrained = isConstrained.sum();
   // coefficient values
   Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> Ak, Bk;
 
-  pvU.resize(numF,2);
-  pvV.resize(numF,2);
-  for (int fi = 0; fi <numF; ++fi)
+  pvU.resize(data.numF,2);
+  pvV.resize(data.numF,2);
+  for (int fi = 0; fi <data.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);
+    const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> &b1 = data.B1.row(fi);
+    const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> &b2 = data.B2.row(fi);
     const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> &u3 = initialSolution.block(fi,0,1,3);
     const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> &v3 = initialSolution.block(fi,3,1,3);
     pvU.row(fi)<< u3.dot(b1), u3.dot(b2);
@@ -537,7 +688,7 @@ solve(const Eigen::VectorXi &isConstrained,
   Ak.resize(numConstrained,1);
   Bk.resize(numConstrained,1);
   int ind = 0;
-  for (int i = 0; i <numF; ++i)
+  for (int i = 0; i <data.numF; ++i)
   {
     if(isConstrained[i])
     {
@@ -546,187 +697,70 @@ 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;
-  
+
   evaluateConjugacy(conjValues);
   meanConj = conjValues.cwiseAbs().mean();
   maxConj = conjValues.cwiseAbs().maxCoeff();
   printf("Initial max non-conjugacy: %.5g\n",maxConj);
-  
-  smoothnessValue = (Acoeff.adjoint()*DDA*Acoeff + Bcoeff.adjoint()*DDB*Bcoeff).real()[0];
+
+  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()*DDA*Acoeff + Bcoeff.adjoint()*DDB*Bcoeff).real()[0];
-    
+
+
+    smoothnessValue = (Acoeff.adjoint()*data.DDA*Acoeff + Bcoeff.adjoint()*data.DDB*Bcoeff).real()[0];
+
     printf("Smoothness: %.5g\n",smoothnessValue);
-    
+
     evaluateConjugacy(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(numF,6);
-  for (int fi=0; fi<numF; ++fi)
+
+  output.setZero(data.numF,6);
+  for (int fi=0; fi<data.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);
+    const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> &b1 = data.B1.row(fi);
+    const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> &b2 = data.B2.row(fi);
     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;
   }
-  
-  
-  
-  return true;
-}
 
+  if (lambdaOut)
+    *lambdaOut = lambda;
 
 
-template<typename DerivedV, typename DerivedF, typename DerivedO>
-IGL_INLINE void igl::ConjugateFFSolver<DerivedV, DerivedF, DerivedO>::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());
-  
-  
+  return true;
 }
 
-template<typename DerivedV, typename DerivedF, typename DerivedO>
-IGL_INLINE void igl::ConjugateFFSolver<DerivedV, DerivedF, DerivedO>::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();
-      
-      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();
-      
-      // 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();
-      
-      // 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();
-      
-      K[eid] = ktemp;
-    }
-  }
-  
-}
 
 
 template <typename DerivedV, typename DerivedF, typename DerivedO>
@@ -735,17 +769,33 @@ IGL_INLINE void igl::conjugate_frame_fields(const Eigen::PlainObjectBase<Derived
                                             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)
+                                            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(V,F);
+  igl::ConjugateFFSolverData<DerivedV, DerivedF> csdata(V, F);
+  igl::ConjugateFFSolver<DerivedV, DerivedF, DerivedO> cs(csdata, maxIter, lambdaOrtho, lambdaInit, lambdaMultFactor, doHardConstraints);
   cs.solve(isConstrained, initialSolution, output);
 }
 
+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::ConjugateFFSolver<DerivedV, DerivedF, DerivedO> cs(csdata, maxIter, lambdaOrtho, lambdaInit, lambdaMultFactor, doHardConstraints);
+  cs.solve(isConstrained, initialSolution, output, lambdaOut);
+}
 
-#ifndef IGL_HEADER_ONLY
+#ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 #endif

include/igl/conjugate_frame_fields.h

@@ -23,6 +23,9 @@ namespace igl {
   // 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,10 +38,22 @@ namespace igl {
                                          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);
+  
 };
 
 
-#ifdef IGL_HEADER_ONLY
+#ifndef IGL_STATIC_LIBRARY
 #include "conjugate_frame_fields.cpp"
 #endif
 

include/igl/edge_topology.cpp

@@ -9,7 +9,7 @@
 #include <algorithm>
 #include "is_edge_manifold.h"
 
-
+template<typename DerivedV, typename DerivedF>
 IGL_INLINE void igl::edge_topology(
                                    const Eigen::PlainObjectBase<DerivedV>& V,
                                    const Eigen::PlainObjectBase<DerivedF>& F,
@@ -33,18 +33,18 @@ IGL_INLINE void igl::edge_topology(
       ETT.push_back(r);
     }
   std::sort(ETT.begin(),ETT.end());
-  
+
   // count the number of edges (assume manifoldness)
   int En = 1; // the last is always counted
   for(unsigned i=0;i<ETT.size()-1;++i)
     if (!((ETT[i][0] == ETT[i+1][0]) && (ETT[i][1] == ETT[i+1][1])))
       ++En;
-  
+
   EV = Eigen::MatrixXi::Constant((int)(En),2,-1);
   FE = Eigen::MatrixXi::Constant((int)(F.rows()),3,-1);
   EF = Eigen::MatrixXi::Constant((int)(En),2,-1);
   En = 0;
-  
+
   for(unsigned i=0;i<ETT.size();++i)
   {
     if (i == ETT.size()-1 ||
@@ -72,10 +72,10 @@ IGL_INLINE void igl::edge_topology(
     }
     ++En;
   }
-  
+
   // Sort the relation EF, accordingly to EV
   // the first one is the face on the left of the edge
-  
+
   for(unsigned i=0; i<EF.rows(); ++i)
   {
     int fid = EF(i,0);
@@ -86,7 +86,7 @@ IGL_INLINE void igl::edge_topology(
       if ((F(fid,j) == EV(i,0)) && (F(fid,(j+1)%3) == EV(i,1)))
         flip = false;
     }
-    
+
     if (flip)
     {
       int tmp = EF(i,0);
@@ -94,7 +94,7 @@ IGL_INLINE void igl::edge_topology(
       EF(i,1) = tmp;
     }
   }
-  
+
 }
 
 #ifdef IGL_STATIC_LIBRARY

include/igl/n_polyvector.cpp

@@ -7,7 +7,7 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 
 #include <igl/n_polyvector.h>
-#include <igl/edgetopology.h>
+#include <igl/edge_topology.h>
 #include <igl/local_basis.h>
 #include <igl/nchoosek.h>
 #include <igl/slice.h>
@@ -77,7 +77,7 @@ numF(_F.rows()),
 n(_n)
 {
 
-  igl::edgetopology(V,F,EV,F2E,E2F);
+  igl::edge_topology(V,F,EV,F2E,E2F);
   numE = EV.rows();
 
 
@@ -492,6 +492,7 @@ IGL_INLINE void igl::n_polyvector(const Eigen::PlainObjectBase<DerivedV> &V,
 }
 
 
-#ifndef IGL_HEADER_ONLY
+#ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
+template void igl::n_polyvector<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> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::Matrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, -1, -1, 0, -1, -1> const&, Eigen::Matrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, -1, -1, 0, -1, -1>&);
 #endif

include/igl/n_polyvector.h

@@ -16,7 +16,7 @@
 namespace igl {
   //todo
   /// 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
@@ -29,11 +29,11 @@ namespace igl {
                                const Eigen::VectorXi &isConstrained,
                                const Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &cfW,
                                Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &output);
-  
+
 };
 
 
-#ifdef IGL_HEADER_ONLY
+#ifndef IGL_STATIC_LIBRARY
 #include "n_polyvector.cpp"
 #endif
 

include/igl/nchoosek.cpp

@@ -53,6 +53,7 @@ IGL_INLINE void igl::nchoosek(int offset,
 
 
 
-#ifndef IGL_HEADER_ONLY
+#ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
+template igl::nchoosek(int, int, int, std::__1::vector<std::__1::vector<int, std::__1::allocator<int> >, std::__1::allocator<std::__1::vector<int, std::__1::allocator<int> > > >&);
 #endif

include/igl/nchoosek.h

@@ -30,7 +30,7 @@ namespace igl {
 }
 
 
-#ifdef IGL_HEADER_ONLY
+#ifndef IGL_STATIC_LIBRARY
 #include "nchoosek.cpp"
 #endif
 

include/igl/polyroots.cpp

@@ -15,22 +15,22 @@ IGL_INLINE void igl::polyRoots(Eigen::Matrix<S, Eigen::Dynamic,1> &polyCoeff, //
 {
   //  degree
   int n = polyCoeff.rows() - 1;
-  
+
   Eigen::Matrix<S, Eigen::Dynamic, 1> d (n,1);
   d = polyCoeff.tail(n)/polyCoeff(0);
-  
+
   Eigen::Matrix<S, Eigen::Dynamic, Eigen::Dynamic> I; I.setIdentity(n-1,n-1);
   Eigen::Matrix<S, Eigen::Dynamic, 1> z; z.setZero(n-1,1);
-  
+
   Eigen::Matrix<S, Eigen::Dynamic, Eigen::Dynamic> a(n,n);
   a<<-d.transpose(),I,z;
   roots = a.eigenvalues();
-  
+
 }
 
 
 
-#ifndef IGL_HEADER_ONLY
+#ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
+template void igl::polyRoots<std::__1::complex<double>, double>(Eigen::Matrix<std::__1::complex<double>, -1, 1, 0, -1, 1>&, Eigen::Matrix<std::__1::complex<double>, -1, 1, 0, -1, 1>&);
 #endif
-

include/igl/polyroots.h

@@ -15,7 +15,7 @@
 namespace igl {
   //todo
   /// 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
@@ -29,7 +29,7 @@ namespace igl {
 }
 
 
-#ifdef IGL_HEADER_ONLY
+#ifndef IGL_STATIC_LIBRARY
 #include "polyroots.cpp"
 #endif
 

tutorial/507_PolyVFField/main.cpp

@@ -29,7 +29,7 @@ bool key_down(igl::Viewer& viewer, unsigned char key, int modifier)
   using namespace std;
   using namespace Eigen;
 
-  if (key <'1' || key >'4')
+  if (key <'0' || key >'4')
     return false;
 
   viewer.clear_mesh();
@@ -39,6 +39,8 @@ bool key_down(igl::Viewer& viewer, unsigned char key, int modifier)
 
   int num = key  - '0';
 
+  if (num == 0)
+    return false;
   // Interpolate
   cerr<<"Interpolating for n = "<<num<<"... ";
   // Interpolated polyVector field
@@ -89,7 +91,7 @@ int main(int argc, char *argv[])
   MatrixXd temp;
   for (int n =0; n<=3; ++n)
   {
-    char cfile[1024]; sprintf(cfile, "/Users/olkido/Desktop/snail%d.dmat",n+1);
+    char cfile[1024]; sprintf(cfile, "../shared/snail%d.dmat",n+1);
 
     igl::readDMAT(cfile,temp);
     if (n == 0)
@@ -102,7 +104,7 @@ int main(int argc, char *argv[])
   igl::Viewer viewer;
 
   // Plot the original mesh with a texture parametrization
-  key_down(viewer,'1',0);
+  key_down(viewer,'0',0);
 
   // Launch the viewer
   viewer.callback_key_down = &key_down;

tutorial/508_ConjugatePVF/CMakeLists.txt

@@ -0,0 +1,11 @@
+cmake_minimum_required(VERSION 2.6)
+project(508_ConjugatePVF)
+
+include("../CMakeLists.shared")
+
+set(SOURCES
+${PROJECT_SOURCE_DIR}/main.cpp
+)
+
+add_executable(${CMAKE_PROJECT_NAME} ${SOURCES} ${SHARED_SOURCES})
+target_link_libraries(${CMAKE_PROJECT_NAME} ${SHARED_LIBRARIES} ${LIBCOMISO_LIBRARY})

tutorial/508_ConjugatePVF/main.cpp

@@ -0,0 +1,213 @@
+#undef IGL_STATIC_LIBRARY
+#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.h>
+#include <igl/conjugate_frame_fields.h>
+#include <stdlib.h>
+#include <igl/readOFF.h>
+#include <igl/jet.h>
+#include <igl/quad_planarity.h>
+#include <igl/planarize_quad_mesh.h>
+
+// Input mesh
+Eigen::MatrixXd V;
+Eigen::MatrixXi F;
+
+// Face barycenters
+Eigen::MatrixXd B;
+
+
+// Quad mesh generated from smooth field
+Eigen::MatrixXd VQS;
+Eigen::MatrixXi FQS;
+Eigen::MatrixXi FQStri;
+
+// Quad mesh generated from conjugate field
+Eigen::MatrixXd VQC;
+Eigen::MatrixXi FQC;
+Eigen::MatrixXi FQCtri;
+Eigen::MatrixXd VQCplan;
+
+// Scale for visualizing the fields
+double global_scale;
+
+// Input constraints
+Eigen::VectorXi isConstrained;
+Eigen::MatrixXd constraints;
+
+Eigen::MatrixXd smooth_pvf;
+Eigen::MatrixXd conjugate_pvf;
+
+igl::ConjugateFFSolverData<Eigen::MatrixXd, Eigen::MatrixXi> *csdata;
+
+int conjIter = 2;
+int totalConjIter = 0;
+double lambdaOrtho = .1;
+double lambdaInit = 100;
+double lambdaMultFactor = 1.01;
+bool doHardConstraints = true;
+
+bool key_down(igl::Viewer& viewer, unsigned char key, int modifier)
+{
+  using namespace std;
+  using namespace Eigen;
+
+  if (key <'1' || key >'6')
+    return false;
+
+  viewer.clear_mesh();
+  viewer.options.show_lines = false;
+  viewer.options.show_texture = false;
+  if (key <= '3')
+  {
+    viewer.set_mesh(V, F);
+    // Highlight in red the constrained faces
+    MatrixXd C = MatrixXd::Constant(F.rows(),3,1);
+    for (unsigned i=0; i<F.rows();++i)
+      if (isConstrained[i])
+        C.row(i) << 1, 0, 0;
+    viewer.set_colors(C);
+  }
+
+  if (key == '1')
+  {
+    // Frame field constraints
+    MatrixXd F1_t = MatrixXd::Zero(F.rows(),3);
+    MatrixXd F2_t = MatrixXd::Zero(F.rows(),3);
+    for (unsigned i=0; i<F.rows();++i)
+      if (isConstrained[i])
+      {
+        F1_t.row(i) = constraints.block(i,0,1,3);
+        F2_t.row(i) = constraints.block(i,3,1,3);
+      }
+    viewer.add_edges (B - global_scale*F1_t, B + global_scale*F1_t , Eigen::RowVector3d(0,0,1));
+    viewer.add_edges (B - global_scale*F2_t, B + global_scale*F2_t , Eigen::RowVector3d(0,0,1));
+
+  }
+  if (key == '2')
+  {
+    viewer.add_edges (B - global_scale*smooth_pvf.block(0,0,F.rows(),3),
+                      B + global_scale*smooth_pvf.block(0,0,F.rows(),3),
+                      Eigen::RowVector3d(0,1,0));
+    viewer.add_edges (B - global_scale*smooth_pvf.block(0,3,F.rows(),3),
+                      B + global_scale*smooth_pvf.block(0,3,F.rows(),3),
+                      Eigen::RowVector3d(0,1,0));
+  }
+
+  if (key == '3')
+  {
+    if (totalConjIter <50)
+    {
+      double lambdaOut;
+      igl::conjugate_frame_fields(*csdata, isConstrained, conjugate_pvf, conjugate_pvf, conjIter, lambdaOrtho, lambdaInit, lambdaMultFactor, doHardConstraints,
+                                  &lambdaOut);
+      totalConjIter += 2;
+      lambdaInit = lambdaOut;
+    }
+    viewer.add_edges (B - global_scale*conjugate_pvf.block(0,0,F.rows(),3),
+                      B + global_scale*conjugate_pvf.block(0,0,F.rows(),3),
+                      Eigen::RowVector3d(0,1,0));
+    viewer.add_edges (B - global_scale*conjugate_pvf.block(0,3,F.rows(),3),
+                      B + global_scale*conjugate_pvf.block(0,3,F.rows(),3),
+                      Eigen::RowVector3d(0,1,0));
+  }
+
+  if (key == '4')
+  {
+    viewer.set_mesh(VQS, FQStri);
+    
+    // show planarity
+    VectorXd planarity;
+    igl::quad_planarity( VQS, FQS, planarity);
+    MatrixXd Ct;
+    igl::jet(planarity, 0, 0.02, Ct);
+    MatrixXd C(FQStri.rows(),3);
+    C << Ct, Ct;
+    viewer.set_colors(C);
+  }
+  
+  if (key == '5')
+  {
+    viewer.set_mesh(VQC, FQCtri);
+    
+    // show planarity
+    VectorXd planarity;
+    igl::quad_planarity( VQC, FQC, planarity);
+    MatrixXd Ct;
+    igl::jet(planarity, 0, 0.02, Ct);
+    MatrixXd C(FQCtri.rows(),3);
+    C << Ct, Ct;
+    viewer.set_colors(C);
+  }
+  
+  if (key == '6')
+  {
+    igl ::planarize_quad_mesh(VQC, FQC, 50, 0.01, VQCplan);
+    viewer.set_mesh(VQCplan, FQCtri);
+
+    // show planarity
+    VectorXd planarity;
+    igl::quad_planarity( VQCplan, FQC, planarity);
+    MatrixXd Ct;
+    igl::jet(planarity, 0, 0.02, Ct);
+    MatrixXd C(FQCtri.rows(),3);
+    C << Ct, Ct;
+    viewer.set_colors(C);
+  }
+
+  return false;
+}
+
+int main(int argc, char *argv[])
+{
+  using namespace Eigen;
+  using namespace std;
+  // Load a mesh in OBJ format
+  igl::readOBJ("../shared/inspired_mesh.obj", V, F);
+
+  // Compute face barycenters
+  igl::barycenter(V, F, B);
+
+  // Compute scale for visualizing fields
+  global_scale =  .2*igl::avg_edge_length(V, F);
+
+  // Load constraints
+  MatrixXd temp;
+  igl::readDMAT("../shared/inspired_mesh.dmat",temp);
+  isConstrained = temp.block(0,0,temp.rows(),1).cast<int>();
+  constraints = temp.block(0,1,temp.rows(),temp.cols()-1);
+
+  // Interpolate to get a smooth field
+  igl::n_polyvector(V, F, isConstrained, constraints, smooth_pvf);
+
+  // Initialize conjugate field with smooth field
+  csdata = new igl::ConjugateFFSolverData<Eigen::MatrixXd,Eigen::MatrixXi>(V,F);
+  conjugate_pvf = smooth_pvf;
+
+  // Load quad mesh generated by smooth field
+  igl::readOFF("../shared/inspired_mesh_quads_Smooth.off", VQS, FQS);
+  FQStri.resize(2*FQS.rows(), 3);
+  FQStri <<  FQS.col(0),FQS.col(1),FQS.col(2),
+             FQS.col(2),FQS.col(3),FQS.col(0);
+
+  // Load quad mesh generated by conjugate field
+  igl::readOFF("../shared/inspired_mesh_quads_Conjugate.off", VQC, FQC);
+  FQCtri.resize(2*FQC.rows(), 3);
+  FQCtri <<  FQC.col(0),FQC.col(1),FQC.col(2),
+             FQC.col(2),FQC.col(3),FQC.col(0);
+
+
+
+  igl::Viewer viewer;
+
+  // Plot the original mesh with a texture parametrization
+  key_down(viewer,'1',0);
+
+  // Launch the viewer
+  viewer.callback_key_down = &key_down;
+  viewer.launch();
+}

tutorial/cmake/FindLIBCOMISO.cmake

@@ -37,7 +37,9 @@ FIND_LIBRARY(LIBCOMISO_LIBRARY NAMES CoMISo
     ${PROJECT_SOURCE_DIR}/../../CoMISo/build/Build/lib/CoMISo/
     ${PROJECT_SOURCE_DIR}/../../../CoMISo/
     ${PROJECT_SOURCE_DIR}/../../../CoMISo/build/Build/lib/CoMISo/
+    /Users/olkido/Documents/igl/MIQ/src/CoMISo/Build
 )
+#message(STATUS "${LIBCOMISO_LIBRARY}")
 
 if(LIBCOMISO_INCLUDE_DIR AND LIBCOMISO_LIBRARY)
 

tutorial/shared/inspired_mesh.dmat.REMOVED.git-id

@@ -0,0 +1 @@
+ef5f214b46f0bf28e206e69f484dd454eac9befe

tutorial/shared/inspired_mesh.obj.REMOVED.git-id

@@ -0,0 +1 @@
+2983d2746f8c29b585ca036ab45c34fb1e638901

tutorial/shared/inspired_mesh_quads_Conjugate.off.REMOVED.git-id

@@ -0,0 +1 @@
+634f371977a3419379f9b673b828f65232ab5c56

tutorial/shared/inspired_mesh_quads_Smooth.off.REMOVED.git-id

@@ -0,0 +1 @@
+8cfb277e82f1e196a90bfd55e68a7f484d0b5275