* new tutorial example for frame fields

* added new frame field functions


Former-commit-id: 9dadbcac64d8c695a016e80b04b3e564b488ca25

include/igl/comiso/frame_field.cpp

@@ -0,0 +1,1203 @@
+#include "frame_field.h"
+
+#include <igl/tt.h>
+#include <igl/edgetopology.h>
+#include <igl/per_face_normals.h>
+#include <igl/comiso/nrosy.h>
+#include <iostream>
+
+namespace igl
+{
+
+class FrameInterpolator
+{
+public:
+  // Constraint type
+  enum ConstraintType { FREE, SOFT, HARD };
+
+  // Init
+  IGL_INLINE FrameInterpolator(const Eigen::MatrixXd& _V, const Eigen::MatrixXi& _F);
+  IGL_INLINE ~FrameInterpolator();
+
+  // Generate the frame field
+  IGL_INLINE void solve(bool hard_const);
+
+  IGL_INLINE void setConstraint(const int fid, const Eigen::VectorXd& v, ConstraintType type = HARD);
+
+  IGL_INLINE void setConstraintProportionalScale(const int fid, const double scale, ConstraintType type = HARD);
+
+  // Set the ratio between smoothness and soft constraints (0 -> smoothness only, 1 -> soft constr only)
+  IGL_INLINE void setSoftAlpha(double alpha);
+
+  // Reset constraints (at least one constraint must be present or solve will fail)
+  IGL_INLINE void resetConstraints();
+
+  // Return the current field
+  IGL_INLINE Eigen::MatrixXd getFieldPerFace();
+
+  // Return the singularities
+  IGL_INLINE Eigen::VectorXd getSingularityIndexPerVertex();
+
+
+  // -------------- This is untested and unstable code
+
+  IGL_INLINE void setConstraint_polar(const int fid, const Eigen::VectorXd& v, ConstraintType type);
+
+  IGL_INLINE void interpolateSymmetric_polar(const bool bilaplacian);
+
+  // Generate the frame field
+  IGL_INLINE void solve_polar(const bool bilaplacian,bool hard_const);
+
+  // Convert the frame field in the canonical representation
+  IGL_INLINE void frame2canonical_polar(const Eigen::MatrixXd& TP, const Eigen::RowVectorXd& v, double& theta, Eigen::VectorXd& S);
+
+  // Convert the canonical representation in a frame field
+  IGL_INLINE void canonical2frame_polar(const Eigen::MatrixXd& TP, const double theta, const Eigen::VectorXd& S, Eigen::RowVectorXd& v);
+
+  IGL_INLINE Eigen::MatrixXd getFieldPerFace_polar();
+
+  IGL_INLINE void PolarDecomposition(Eigen::MatrixXd V, Eigen::MatrixXd& U, Eigen::MatrixXd& P);
+
+  // Symmetric
+  Eigen::MatrixXd S_polar;
+  std::vector<ConstraintType> S_polar_c;
+
+  // -------------------------------------------------
+
+  // Face Topology
+  Eigen::MatrixXi TT, TTi;
+
+  // Two faces are consistent if their representative vector are taken modulo PI
+  std::vector<bool> edge_consistency;
+  Eigen::MatrixXi   edge_consistency_TT;
+
+private:
+  IGL_INLINE double mod2pi(double d);
+  IGL_INLINE double modpi2(double d);
+  IGL_INLINE double modpi(double d);
+
+  // Convert a direction on the tangent space into an angle
+  IGL_INLINE double vector2theta(const Eigen::MatrixXd& TP, const Eigen::RowVectorXd& v);
+
+  // Convert an angle in a vector in the tangent space
+  IGL_INLINE Eigen::RowVectorXd theta2vector(const Eigen::MatrixXd& TP, const double theta);
+
+  // Convert the frame field in the canonical representation
+  IGL_INLINE void frame2canonical(const Eigen::MatrixXd& TP, const Eigen::RowVectorXd& v, double& theta, double& beta, Eigen::RowVectorXd& d);
+
+  // Convert the canonical representation in a frame field
+  IGL_INLINE void canonical2frame(const Eigen::MatrixXd& TP, const double theta, const double beta, const Eigen::RowVectorXd& d, Eigen::RowVectorXd& v);
+
+  // Interpolate the cross field (theta)
+  IGL_INLINE void interpolateCross(bool hard_const);
+
+  // Interpolate the skewness (beta)
+  IGL_INLINE void interpolateSkewness();
+
+  // Interpolate the scale (d)
+  IGL_INLINE void interpolateScale();
+
+  // Compute difference between reference frames
+  IGL_INLINE void computek();
+
+  // Compute edge consistency
+  IGL_INLINE void compute_edge_consistency();
+
+  // Cross field direction
+  Eigen::VectorXd thetas;
+  std::vector<ConstraintType> thetas_c;
+
+  // Skewness
+  Eigen::VectorXd betas;
+  std::vector<ConstraintType> betas_c;
+
+  // Scale
+  Eigen::MatrixXd ds;
+  std::vector<ConstraintType> ds_c;
+
+  // Soft constraints strenght
+  double softAlpha;
+
+  // Edge Topology
+  Eigen::MatrixXi EV, FE, EF;
+  std::vector<bool> isBorderEdge;
+
+  // Angle between two reference frames
+  // R(k) * t0 = t1
+  Eigen::VectorXd k;
+
+  // Mesh
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+
+  // Normals per face
+  Eigen::MatrixXd N;
+
+  // Singularity index
+  Eigen::VectorXd singularityIndex;
+
+  // Reference frame per triangle
+  std::vector<Eigen::MatrixXd> TPs;
+
+};
+
+
+
+
+FrameInterpolator::FrameInterpolator(const Eigen::MatrixXd& _V, const Eigen::MatrixXi& _F)
+{
+  using namespace std;
+  using namespace Eigen;
+
+  V = _V;
+  F = _F;
+
+  assert(V.rows() > 0);
+  assert(F.rows() > 0);
+
+
+  // Generate topological relations
+  igl::tt(V,F,TT,TTi);
+  igl::edgetopology(V,F, EV, FE, EF);
+
+  // Flag border edges
+  isBorderEdge.resize(EV.rows());
+  for(unsigned i=0; i<EV.rows(); ++i)
+    isBorderEdge[i] = (EF(i,0) == -1) || ((EF(i,1) == -1));
+
+  // Generate normals per face
+  igl::per_face_normals(V, F, N);
+
+  // Generate reference frames
+  for(unsigned fid=0; fid<F.rows(); ++fid)
+  {
+    // First edge
+    Vector3d e1 = V.row(F(fid,1)) - V.row(F(fid,0));
+    e1.normalize();
+    Vector3d e2 = N.row(fid);
+    e2 = e2.cross(e1);
+    e2.normalize();
+
+    MatrixXd TP(2,3);
+    TP << e1.transpose(), e2.transpose();
+    TPs.push_back(TP);
+  }
+
+  // Reset the constraints
+  resetConstraints();
+
+  // Compute k, differences between reference frames
+  computek();
+
+  softAlpha = 0.5;
+
+  // Alloc internal variables
+  thetas            = VectorXd::Zero(F.rows());
+//  thetas_c.resize(F.rows());
+
+  betas            = VectorXd::Zero(F.rows());
+//  betas_c.resize(F.rows());
+
+  ds                = MatrixXd::Constant(F.rows(),2,1);
+//  ds_c.resize(F.rows());
+
+  S_polar = MatrixXd::Zero(F.rows(),3);
+//  S_polar_c.resize(F.rows());
+
+  singularityIndex = VectorXd::Zero(V.rows());
+
+  compute_edge_consistency();
+
+}
+
+FrameInterpolator::~FrameInterpolator()
+{
+
+}
+
+void FrameInterpolator::solve(bool hard_const)
+{
+    interpolateCross(hard_const);
+    interpolateSkewness();
+    interpolateScale();
+}
+
+void FrameInterpolator::setConstraint(const int fid, const Eigen::VectorXd& v, ConstraintType type)
+{
+  using namespace std;
+  using namespace Eigen;
+  
+  double   t_;
+  double   b_;
+  RowVectorXd d_;
+  frame2canonical(TPs[fid],v,t_,b_,d_);
+
+//  cerr << "TP: " << endl << TPs[fid] << endl;
+//  cerr << "v: " << endl << TPs[fid] << endl;
+//  cerr << "t_: " << endl << TPs[fid] << endl;
+//  cerr << "b_: " << endl << TPs[fid] << endl;
+//  cerr << "d_: " << endl << TPs[fid] << endl;
+
+  Eigen::RowVectorXd v2;
+  canonical2frame(TPs[fid], t_, b_, d_, v2);
+//  cerr << "Cosntraint: " << t_ << " " << b_ << " " << d_ << endl;
+//  cerr << "Compare: " << endl << v.transpose() << endl << v2 << endl;
+
+  thetas(fid)   = t_;
+  thetas_c[fid] = type;
+
+  betas(fid)    = b_;
+  betas_c[fid]  = type;
+
+  ds.row(fid)   = d_;
+  ds_c[fid]     = type;
+}
+
+void FrameInterpolator::setConstraintProportionalScale(const int fid, const double scale, ConstraintType type)
+{
+  using namespace std;
+  using namespace Eigen;
+  
+
+  assert(scale>0);
+  ds.row(fid)   = ds.row(fid) * scale;
+  ds_c[fid]     = type;
+}
+
+void FrameInterpolator::setSoftAlpha(double alpha)
+{
+  assert(alpha >= 0 && alpha <= 1);
+  softAlpha = alpha;
+}
+
+
+Eigen::MatrixXd FrameInterpolator::getFieldPerFace()
+{
+  Eigen::MatrixXd R(F.rows(),6);
+  for (unsigned i=0; i<F.rows(); ++i)
+  {
+    Eigen::RowVectorXd v;
+    canonical2frame(TPs[i],thetas(i),betas(i),ds.row(i),v);
+    R.row(i) = v;
+  }
+  return R;
+}
+
+Eigen::VectorXd FrameInterpolator::getSingularityIndexPerVertex()
+{
+  assert(0);
+}
+
+
+double FrameInterpolator::mod2pi(double d)
+{
+  while(d<0)
+    d = d + (2.0*M_PI);
+
+  return fmod(d, (2.0*M_PI));
+}
+
+double FrameInterpolator::modpi2(double d)
+{
+  while(d<0)
+    d = d + (M_PI/2.0);
+
+  return fmod(d, (M_PI/2.0));
+}
+
+double FrameInterpolator::modpi(double d)
+{
+  while(d<0)
+    d = d + (M_PI);
+
+  return fmod(d, (M_PI));
+}
+
+
+double FrameInterpolator::vector2theta(const Eigen::MatrixXd& TP, const Eigen::RowVectorXd& v)
+{
+  // Project onto the tangent plane
+  Eigen::Vector2d vp = TP * v.transpose();
+
+  // Convert to angle
+  double theta = atan2(vp(1),vp(0));
+
+//  cerr << v << endl << theta2vector(TP, theta) << endl;
+//  assert((v.normalized() - theta2vector(TP, theta)).norm() < 10e-5);
+
+  return theta;
+}
+
+Eigen::RowVectorXd FrameInterpolator::theta2vector(const Eigen::MatrixXd& TP, const double theta)
+{
+  Eigen::Vector2d vp(cos(theta),sin(theta));
+  return vp.transpose() * TP;
+}
+
+void FrameInterpolator::frame2canonical(const Eigen::MatrixXd& TP, const Eigen::RowVectorXd& v, double& theta, double& beta, Eigen::RowVectorXd& d)
+{
+  using namespace std;
+  using namespace Eigen;
+  
+
+  bool debug = false;
+  if (debug)
+  {
+    cerr << "TP:" << endl << TP << endl;
+    cerr << "v:" << endl << v << endl;
+  }
+
+  RowVectorXd v0 = v.segment<3>(0);
+  RowVectorXd v1 = v.segment<3>(3);
+
+  if (debug)
+  {
+    cerr << "v0:" << endl << v0 << endl;
+    cerr << "v1:" << endl << v1 << endl;
+  }
+
+  // Find canonical cross field
+  double theta0 = mod2pi(vector2theta(TP,v0));
+  double theta1 = mod2pi(vector2theta(TP,v1));
+  //cerr << mod2pi(theta0) << " " << mod2pi(theta1) << endl;
+
+  // Find the closest pair
+//  if (mod2pi(abs(theta0 + M_PI - theta1)) < mod2pi(abs(theta0 - theta1)))
+//      theta0 = theta0 + M_PI;
+
+  // Express theta1 as a displacement wrt theta0
+  theta1 = mod2pi(theta1 - theta0);
+
+  // Find the closest representative
+  if (theta1 > M_PI)
+    theta1 -= M_PI;
+  assert(theta1 < M_PI);
+
+  // Extract theta and beta
+  theta = mod2pi(theta1/2.0 + theta0);
+  beta  = theta1/2.0;
+
+  //cerr << mod2pi(beta + theta) << " " << mod2pi(-beta + theta) << endl;
+
+  d.resize(2);
+  // Find the scaling factors
+  d(0) = v0.norm();
+  d(1) = v1.norm();
+
+  if (debug)
+  {
+    cerr << "d:" << endl << d << endl;
+    cerr << "thetavector:" << endl << theta2vector(TP, theta) << endl;
+  }
+
+}
+
+void FrameInterpolator::canonical2frame(const Eigen::MatrixXd& TP, const double theta, const double beta, const Eigen::RowVectorXd& d, Eigen::RowVectorXd& v)
+{
+  using namespace std;
+  using namespace Eigen;
+  
+
+  assert(d.cols() == 2 && d.rows() == 1);
+
+  double theta0 = theta + M_PI - beta;
+  double theta1 = theta + M_PI + beta;
+
+  RowVectorXd v0 = theta2vector(TP,modpi(theta0)) * d(0);
+  RowVectorXd v1 = theta2vector(TP,modpi(theta1)) * d(1);
+
+//  cerr << "v0 : " << v0 << endl;
+//  cerr << "v1 : " << v1 << endl;
+
+  v.resize(6);
+  v << v0, v1;
+}
+
+void FrameInterpolator::interpolateCross(bool hard_const)
+{
+  using namespace std;
+  using namespace Eigen;
+
+  NRosyField nrosy(V,F);
+  nrosy.setSoftAlpha(softAlpha);
+
+  for (unsigned i=0; i<F.rows(); ++i)
+  {
+    switch (thetas_c[i])
+    {
+      case FREE:
+      break;
+      case SOFT:
+//      nrosy.setConstraintSoft(i,1,theta2vector(TPs[i],thetas(i)));
+//      cerr << theta2vector(TPs[i],thetas(i)) << endl;
+//      break;
+      case HARD:
+//      nrosy.setConstraintHard(i,theta2vector(TPs[i],thetas(i)));
+        if (hard_const)
+          nrosy.setConstraintHard(i,theta2vector(TPs[i],thetas(i)));
+        else
+          nrosy.setConstraintSoft(i,1,theta2vector(TPs[i],thetas(i)));
+      break;
+      default:
+      assert(0);
+    }
+  }
+
+  nrosy.solve(4);
+
+  MatrixXd R = nrosy.getFieldPerFace();
+//  std::cerr << "Cross: " << R << std::endl;
+
+  assert(R.rows() == F.rows());
+  for (unsigned i=0; i<F.rows(); ++i)
+    if (thetas_c[i] != HARD)
+      thetas(i) = vector2theta(TPs[i],R.row(i));
+}
+
+void FrameInterpolator::interpolateSkewness()
+{
+  using namespace std;
+  using namespace Eigen;
+  
+  compute_edge_consistency();
+
+  // Generate enriched Laplacian matrix (see notes)
+  typedef Eigen::Triplet<double> triplet;
+  std::vector<triplet> triplets;
+  triplets.reserve(4*F.rows());
+
+  VectorXd b = VectorXd::Zero(F.rows());
+
+  // Build L and b
+  for (unsigned eid=0; eid<EF.rows(); ++eid)
+  {
+    if (!isBorderEdge[eid])
+    {
+      int i = EF(eid,0);
+      int j = EF(eid,1);
+      double v = edge_consistency[eid]? 1 : -1;
+
+      // Off-diagonal, symmetric
+      triplets.push_back(triplet(i,j,v));
+      triplets.push_back(triplet(j,i,v));
+      // Diagonal
+      triplets.push_back(triplet(i,i,-1));
+      triplets.push_back(triplet(j,j,-1));
+
+      if (!edge_consistency[eid])
+      {
+        b(i) -= M_PI/2.0; // CHECK
+        b(j) -= M_PI/2.0; // CHECK
+      }
+    }
+  }
+
+  // Add soft constraints
+  double w = 100000;
+  for (unsigned fid=0; fid < F.rows(); ++fid)
+  {
+    if (betas_c[fid] != FREE)
+    {
+      triplets.push_back(triplet(fid,fid,w));
+      b(fid) += w*betas(fid);
+    }
+  }
+
+  SparseMatrix<double> L(F.rows(),F.rows());
+  L.setFromTriplets(triplets.begin(), triplets.end());
+
+  // Solve Lx = b;
+
+  SimplicialLDLT<SparseMatrix<double> > solver;
+  solver.compute(L);
+
+  if(solver.info()!=Success)
+  {
+    std::cerr << "Cholesky failed - frame_interpolator.cpp" << std::endl;
+    assert(0);
+  }
+
+  VectorXd x;
+  x = solver.solve(b);
+
+//  cerr << "Skewness: " << x << endl;
+
+  if(solver.info()!=Success)
+  {
+    std::cerr << "Linear solve failed - frame_interpolator.cpp" << std::endl;
+    assert(0);
+  }
+
+  // Copy back the result
+  betas = x;
+}
+
+void FrameInterpolator::interpolateScale()
+{
+  using namespace std;
+  using namespace Eigen;
+
+  compute_edge_consistency();
+
+  // Generate enriched Laplacian matrix (see notes)
+  // the variables here are d1, d2
+  typedef Eigen::Triplet<double> triplet;
+  std::vector<triplet> triplets;
+  triplets.reserve(4*F.rows()*2);
+
+  VectorXd b = VectorXd::Zero(F.rows()*2);
+
+  // Build L and b
+  for (unsigned eid=0; eid<EF.rows(); ++eid)
+  {
+    if (!isBorderEdge[eid])
+    {
+      int d1  = EF(eid,0);
+      int d2  = EF(eid,0)+F.rows();
+      int d1p = EF(eid,1);
+      int d2p = EF(eid,1)+F.rows();
+
+      if (edge_consistency[eid])
+      {
+        triplets.push_back(triplet(d1 ,d1p,1));
+        triplets.push_back(triplet(d2 ,d2p,1));
+        triplets.push_back(triplet(d1p,d1 ,1));
+        triplets.push_back(triplet(d2p,d2 ,1));
+      }
+      else
+      {
+        triplets.push_back(triplet(d1 ,d2p,1));
+        triplets.push_back(triplet(d2 ,d1p,1));
+        triplets.push_back(triplet(d1p,d2 ,1));
+        triplets.push_back(triplet(d2p,d1 ,1));
+      }
+
+      // Diagonal
+      triplets.push_back(triplet(d1,d1,  -1));
+      triplets.push_back(triplet(d2,d2,  -1));
+      triplets.push_back(triplet(d1p,d1p,-1));
+      triplets.push_back(triplet(d2p,d2p,-1));
+    }
+  }
+
+  // Add soft constraints
+  double w = 100000;
+  for (unsigned fid=0; fid < F.rows(); ++fid)
+  {
+    if (ds_c[fid] != FREE)
+    {
+      int d1 = fid;
+      int d2 = fid + F.rows();
+
+      triplets.push_back(triplet(d1,d1,w));
+      b(d1) += w*ds(fid,0);
+
+      triplets.push_back(triplet(d2,d2,w));
+      b(d2) += w*ds(fid,1);
+    }
+  }
+
+  SparseMatrix<double> L(F.rows()*2,F.rows()*2);
+  L.setFromTriplets(triplets.begin(), triplets.end());
+
+  // Solve Lx = b;
+
+  SimplicialLDLT<SparseMatrix<double> > solver;
+  solver.compute(L);
+
+  if(solver.info()!=Success)
+  {
+    std::cerr << "Cholesky failed - frame_interpolator.cpp" << std::endl;
+    assert(0);
+  }
+
+  VectorXd x;
+  x = solver.solve(b);
+
+//  cerr << "Scale: " << endl << x << endl;
+  if(solver.info()!=Success)
+  {
+    std::cerr << "Linear solve failed - frame_interpolator.cpp" << std::endl;
+    assert(0);
+  }
+
+  // Copy back the result
+  ds << x.block(0, 0, ds.rows(), 1), x.block(ds.rows(), 0, ds.rows(), 1);
+}
+
+void FrameInterpolator::resetConstraints()
+{
+  thetas_c.resize(F.rows());
+  betas_c.resize(F.rows());
+  ds_c.resize(F.rows());
+  S_polar_c.resize(F.rows());
+
+  for(unsigned i=0; i<F.rows(); ++i)
+  {
+    thetas_c[i]  = FREE;
+    betas_c[i]   = FREE;
+    ds_c[i]      = FREE;
+    S_polar_c[i] = FREE;
+  }
+
+}
+
+void FrameInterpolator::compute_edge_consistency()
+{
+  using namespace std;
+  using namespace Eigen;
+
+  // Compute per-edge consistency
+  edge_consistency.resize(EF.rows());
+  edge_consistency_TT = MatrixXi::Constant(TT.rows(),3,-1);
+
+  // For every non-border edge
+  for (unsigned eid=0; eid<EF.rows(); ++eid)
+  {
+    if (!isBorderEdge[eid])
+    {
+      int fid0 = EF(eid,0);
+      int fid1 = EF(eid,1);
+
+      double theta0 = thetas(fid0);
+      double theta1 = thetas(fid1);
+
+      theta0 = theta0 + k(eid);
+
+      double r = modpi(theta0-theta1);
+
+      edge_consistency[eid] = r < M_PI/4.0 || r > 3*(M_PI/4.0);
+
+      // Copy it into edge_consistency_TT
+      int i1 = -1;
+      int i2 = -1;
+      for (unsigned i=0; i<3; ++i)
+      {
+        if (TT(fid0,i) == fid1)
+          i1 = i;
+        if (TT(fid1,i) == fid0)
+          i2 = i;
+      }
+      assert(i1 != -1);
+      assert(i2 != -1);
+
+      edge_consistency_TT(fid0,i1) = edge_consistency[eid];
+      edge_consistency_TT(fid1,i2) = edge_consistency[eid];
+    }
+  }
+}
+
+void FrameInterpolator::computek()
+{
+  using namespace std;
+  using namespace Eigen;
+
+  k.resize(EF.rows());
+
+  // For every non-border edge
+  for (unsigned eid=0; eid<EF.rows(); ++eid)
+  {
+    if (!isBorderEdge[eid])
+    {
+      int fid0 = EF(eid,0);
+      int fid1 = EF(eid,1);
+
+      Vector3d N0 = N.row(fid0);
+      //Vector3d N1 = N.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 (EV(eid,0) == F(fid0,i))
+          fid0_vc = i;
+        if (EV(eid,1) == F(fid1,i))
+          fid1_vc = i;
+      }
+      assert(fid0_vc != -1);
+      assert(fid1_vc != -1);
+
+      Vector3d 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
+      MatrixXd P(3,3);
+      VectorXd o = V.row(F(fid0,fid0_vc));
+      VectorXd tmp = -N0.cross(common_edge);
+      P << common_edge, tmp, N0;
+      P.transposeInPlace();
+
+
+      MatrixXd V0(3,3);
+      V0.row(0) = V.row(F(fid0,0)).transpose() -o;
+      V0.row(1) = V.row(F(fid0,1)).transpose() -o;
+      V0.row(2) = V.row(F(fid0,2)).transpose() -o;
+
+      V0 = (P*V0.transpose()).transpose();
+
+      assert(V0(0,2) < 10e-10);
+      assert(V0(1,2) < 10e-10);
+      assert(V0(2,2) < 10e-10);
+
+      MatrixXd V1(3,3);
+      V1.row(0) = V.row(F(fid1,0)).transpose() -o;
+      V1.row(1) = V.row(F(fid1,1)).transpose() -o;
+      V1.row(2) = V.row(F(fid1,2)).transpose() -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));
+
+      MatrixXd R(3,3);
+      R << 1,          0,            0,
+           0, cos(alpha), -sin(alpha) ,
+           0, sin(alpha),  cos(alpha);
+      V1 = (R*V1.transpose()).transpose();
+
+      assert(V1(0,2) < 10e-10);
+      assert(V1(1,2) < 10e-10);
+      assert(V1(2,2) < 10e-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
+      VectorXd ref0 = V0.row(1) - V0.row(0);
+      VectorXd 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...
+      MatrixXd R2(2,2);
+      R2 << cos(ktemp), -sin(ktemp), sin(ktemp), cos(ktemp);
+
+      tmp = R2*ref0.head<2>();
+
+      assert(tmp(0) - ref1(0) < (0.000001));
+      assert(tmp(1) - ref1(1) < (0.000001));
+
+      k[eid] = ktemp;
+    }
+  }
+
+}
+
+
+  void FrameInterpolator::frame2canonical_polar(const Eigen::MatrixXd& TP, const Eigen::RowVectorXd& v, double& theta, Eigen::VectorXd& S_v)
+{
+  using namespace std;
+  using namespace Eigen;
+
+  RowVectorXd v0 = v.segment<3>(0);
+  RowVectorXd v1 = v.segment<3>(3);
+
+  // Project onto the tangent plane
+  Vector2d vp0 = TP * v0.transpose();
+  Vector2d vp1 = TP * v1.transpose();
+
+  // Assemble matrix
+  MatrixXd M(2,2);
+  M << vp0, vp1;
+
+  if (M.determinant() < 0)
+    M.col(1) = -M.col(1);
+
+  assert(M.determinant() > 0);
+
+  // cerr << "M: " << M << endl;
+
+  MatrixXd R,S;
+  // PolarDecomposition
+  PolarDecomposition(M,R,S);
+
+//  cerr << M << endl;
+//  cerr << "--------" << endl;
+//  cerr << R*S << endl;
+
+  // Finally, express the cross field as an angle
+  theta = atan2(R(1,0),R(0,0)); // atan2(R(1,0),R(0,0));
+
+  MatrixXd R2(2,2);
+  R2 << cos(theta), -sin(theta), sin(theta), cos(theta);
+
+//  cerr << "R: " << R << endl;
+//  cerr << "Should be zero2: " << R2-R << endl;
+
+  assert((R2-R).norm() < 10e-8);
+
+
+  // Convert into rotation invariant form
+  S = R * S * R.inverse();
+
+  // cerr << "R" << endl << R << endl;
+  // cerr << "S" << endl << S << endl;
+
+  // Copy in vector form
+  S_v = VectorXd(3);
+  S_v << S(0,0), S(0,1), S(1,1);
+
+}
+
+  void FrameInterpolator::canonical2frame_polar(const Eigen::MatrixXd& TP, const double theta, const Eigen::VectorXd& S_v, Eigen::RowVectorXd& v)
+{
+  using namespace std;
+  using namespace Eigen;
+
+  assert(S_v.size() == 3);
+
+  MatrixXd S(2,2);
+  S << S_v(0), S_v(1), S_v(1), S_v(2);
+
+  // Convert angle in vector in the tangent plane
+  // Vector2d vp(cos(theta),sin(theta));
+
+  // First reconstruct R
+  MatrixXd R(2,2);
+  // R.col(0) << cos(theta), sin(theta);
+  // R(0,1) = -R(1,0);
+  // R(1,1) =  R(0,0);
+  R << cos(theta), -sin(theta), sin(theta), cos(theta);
+
+
+  // cerr << "R2" << endl << R << endl;
+  // cerr << "S2" << endl << S << endl;
+
+  // Multiply to reconstruct
+  //MatrixXd M = R * S;
+
+  // Rotation invariant reconstruction
+  MatrixXd M = S * R;
+
+  // cerr << "M2: " << M << endl;
+
+  Vector2d vp0(M(0,0),M(1,0));
+  Vector2d vp1(M(0,1),M(1,1));
+
+  // Unproject the vectors
+  RowVectorXd v0 = vp0.transpose() * TP;
+  RowVectorXd v1 = vp1.transpose() * TP;
+
+  v.resize(6);
+  v << v0, v1;
+}
+
+void FrameInterpolator::solve_polar(const bool bilaplacian, bool hard_const)
+{
+    interpolateCross(hard_const);
+    //interpolateSymmetric_polar(bilaplacian);
+  interpolateSymmetric_polar(false);
+}
+
+void FrameInterpolator::interpolateSymmetric_polar(const bool bilaplacian)
+{
+  using namespace std;
+  using namespace Eigen;
+
+  //compute_edge_consistency();
+
+  // Generate uniform Laplacian matrix (see notes)
+  typedef Eigen::Triplet<double> triplet;
+  std::vector<triplet> triplets;
+
+  // Variables are stacked as x1,y1,z1,x2,y2,z2
+  triplets.reserve(3*4*F.rows());
+
+  MatrixXd b = MatrixXd::Zero(3*F.rows(),1);
+
+  // Build L and b
+  for (unsigned eid=0; eid<EF.rows(); ++eid)
+  {
+    if (!isBorderEdge[eid])
+    {
+      for (int z=0;z<2;++z)
+      {
+        // W = [w_a, w_b
+        //      w_b, w_c]
+        //
+
+        // It is not symmetric
+        int i    = EF(eid,z==0?0:1);
+        int j    = EF(eid,z==0?1:0);
+
+        int w_a_0 = (i*3)+0;
+        int w_b_0 = (i*3)+1;
+        int w_c_0 = (i*3)+2;
+
+        int w_a_1 = (j*3)+0;
+        int w_b_1 = (j*3)+1;
+        int w_c_1 = (j*3)+2;
+
+        // Rotation to change frame
+        double r_a =  cos(z==1?k(eid):-k(eid));
+        double r_b = -sin(z==1?k(eid):-k(eid));
+        double r_c =  sin(z==1?k(eid):-k(eid));
+        double r_d =  cos(z==1?k(eid):-k(eid));
+
+        if (true)
+        {
+          if (true)
+          {
+            // First term
+            // w_a_0 = r_a^2 w_a_1 + 2 r_a r_b w_b_1 + r_b^2 w_c_1 = 0
+            triplets.push_back(triplet(w_a_0,w_a_0,                -1 ));
+            triplets.push_back(triplet(w_a_0,w_a_1,           r_a*r_a ));
+            triplets.push_back(triplet(w_a_0,w_b_1,       2 * r_a*r_b ));
+            triplets.push_back(triplet(w_a_0,w_c_1,           r_b*r_b ));
+
+            // Second term
+            // w_b_0 = r_a r_c w_a + (r_b r_c + r_a r_d) w_b + r_b r_d w_c
+            triplets.push_back(triplet(w_b_0,w_b_0,                -1 ));
+            triplets.push_back(triplet(w_b_0,w_a_1,           r_a*r_c ));
+            triplets.push_back(triplet(w_b_0,w_b_1, r_b*r_c + r_a*r_d ));
+            triplets.push_back(triplet(w_b_0,w_c_1,           r_b*r_d ));
+
+            // Third term
+            // w_c_0 = r_c^2 w_a + 2 r_c r_d w_b +  r_d^2 w_c
+            triplets.push_back(triplet(w_c_0,w_c_0,                -1 ));
+            triplets.push_back(triplet(w_c_0,w_a_1,           r_c*r_c ));
+            triplets.push_back(triplet(w_c_0,w_b_1,       2 * r_c*r_d ));
+            triplets.push_back(triplet(w_c_0,w_c_1,           r_d*r_d ));
+          }
+          else
+          {
+            // First term
+            // w_a_0 = r_a^2 w_a_1 + 2 r_a r_b w_b_1 + r_b^2 w_c_1 = 0
+            triplets.push_back(triplet(w_a_0,w_a_1,           r_a*r_a ));
+            triplets.push_back(triplet(w_a_0,w_b_1,       2 * r_a*r_b ));
+            triplets.push_back(triplet(w_a_0,w_c_1,           r_b*r_b ));
+            triplets.push_back(triplet(w_a_0,w_a_0, -          r_a*r_a ));
+            triplets.push_back(triplet(w_a_0,w_a_0, -      2 * r_a*r_b ));
+            triplets.push_back(triplet(w_a_0,w_a_0, -          r_b*r_b ));
+
+            // Second term
+            // w_b_0 = r_a r_c w_a + (r_b r_c + r_a r_d) w_b + r_b r_d w_c
+            triplets.push_back(triplet(w_b_0,w_a_1,           r_a*r_c ));
+            triplets.push_back(triplet(w_b_0,w_b_1, r_b*r_c + r_a*r_d ));
+            triplets.push_back(triplet(w_b_0,w_c_1,           r_b*r_d ));
+            triplets.push_back(triplet(w_b_0,w_b_0, -          r_a*r_c ));
+            triplets.push_back(triplet(w_b_0,w_b_0, -r_b*r_c + r_a*r_d ));
+            triplets.push_back(triplet(w_b_0,w_b_0, -          r_b*r_d ));
+
+            // Third term
+            // w_c_0 = r_c^2 w_a + 2 r_c r_d w_b +  r_d^2 w_c
+            triplets.push_back(triplet(w_c_0,w_a_1,           r_c*r_c ));
+            triplets.push_back(triplet(w_c_0,w_b_1,       2 * r_c*r_d ));
+            triplets.push_back(triplet(w_c_0,w_c_1,           r_d*r_d ));
+            triplets.push_back(triplet(w_c_0,w_c_0, -          r_c*r_c ));
+            triplets.push_back(triplet(w_c_0,w_c_0, -      2 * r_c*r_d ));
+            triplets.push_back(triplet(w_c_0,w_c_0, -          r_d*r_d ));
+
+          }
+        }
+        else
+        {
+          // Simple Laplacian
+          // w_a_0 = r_a^2 w_a_1 + 2 r_a r_b w_b_1 + r_b^2 w_c_1 = 0
+          triplets.push_back(triplet(w_a_0,w_a_0,                -1 ));
+          triplets.push_back(triplet(w_a_0,w_a_1,                 1 ));
+
+          // Second term
+          // w_b_0 = r_a r_c w_a + (r_b r_c + r_a r_d) w_b + r_b r_d w_c
+          triplets.push_back(triplet(w_b_0,w_b_0,                -1 ));
+          triplets.push_back(triplet(w_b_0,w_b_1,                 1 ));
+
+          // Third term
+          // w_c_0 = r_c^2 w_a + 2 r_c r_d w_b +  r_d^2 w_c
+          triplets.push_back(triplet(w_c_0,w_c_0,                -1 ));
+          triplets.push_back(triplet(w_c_0,w_c_1,                 1 ));
+        }
+      }
+    }
+  }
+
+  SparseMatrix<double> L(3*F.rows(),3*F.rows());
+  L.setFromTriplets(triplets.begin(), triplets.end());
+
+  triplets.clear();
+
+  // Add soft constraints
+  double w = 100000;
+  for (unsigned fid=0; fid < F.rows(); ++fid)
+  {
+    if (S_polar_c[fid] != FREE)
+    {
+      for (unsigned i=0;i<3;++i)
+      {
+        triplets.push_back(triplet(3*fid + i,3*fid + i,w));
+        b(3*fid + i) += w*S_polar(fid,i);
+      }
+    }
+  }
+
+  SparseMatrix<double> soft(3*F.rows(),3*F.rows());
+  soft.setFromTriplets(triplets.begin(), triplets.end());
+
+  SparseMatrix<double> M;
+
+  if (!bilaplacian)
+    M = L + soft;
+  else
+    M = L * L + soft;
+
+  // Solve Lx = b;
+
+//  MatrixXd Dense(L);
+//
+//  cerr << "Is sym:" << (Dense - Dense.transpose()).maxCoeff() << endl;
+
+//  SimplicialLDLT<SparseMatrix<double> > solver;
+  SparseLU<SparseMatrix<double> > solver;
+
+  solver.compute(M);
+
+  if(solver.info()!=Success)
+  {
+    std::cerr << "Cholesky failed - frame_interpolator.cpp" << std::endl;
+    assert(0);
+  }
+
+  MatrixXd x;
+  x = solver.solve(b);
+
+//  cerr << "Skewness: " << x << endl;
+
+  if(solver.info()!=Success)
+  {
+    std::cerr << "Linear solve failed - frame_interpolator.cpp" << std::endl;
+    assert(0);
+  }
+
+  S_polar = MatrixXd::Zero(F.rows(),3);
+
+  // Copy back the result
+  for (unsigned i=0;i<F.rows();++i)
+    S_polar.row(i) << x(i*3+0), x(i*3+1), x(i*3+2);
+
+  for (unsigned i=0;i<F.rows();++i)
+  {
+    if (false)
+    {
+      MatrixXd T_POLAR(2,2);
+      T_POLAR << x(i*3+0), x(i*3+1), x(i*3+1), x(i*3+2);
+
+      Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> S_es(T_POLAR);
+      VectorXd eivals = S_es.eigenvalues();
+
+//      cerr << "eivals:" << eivals << endl;
+      double t = eivals.minCoeff();
+      if (t < 0)
+      {
+        cerr << "NOOOOOOOO!" << endl;
+        exit(1);
+      }
+//      cerr << S_polar << endl;
+    }
+  }
+}
+
+void FrameInterpolator::setConstraint_polar(const int fid, const Eigen::VectorXd& v, ConstraintType type)
+{
+  using namespace std;
+  using namespace Eigen;
+
+  double   t_;
+  VectorXd S_;
+
+  frame2canonical_polar(TPs[fid],v,t_,S_);
+
+  Eigen::RowVectorXd v2;
+  canonical2frame_polar(TPs[fid], t_, S_, v2);
+
+
+//  cerr << "Constraint: " << t_ << " " << S_ << endl;
+//  cerr << "Compare:    " << endl << v.transpose() << endl << v2 << endl;
+//  cerr << "Diff:    " << endl << v.transpose() - v2 << endl;
+//  cerr << "Diff:    " << endl << v.transpose() + v2 << endl;
+
+  thetas(fid)   = t_;
+  thetas_c[fid] = type;
+
+  S_polar.row(fid) = S_;
+  S_polar_c[fid]   = type;
+
+}
+
+Eigen::MatrixXd FrameInterpolator::getFieldPerFace_polar()
+{
+  using namespace std;
+  using namespace Eigen;
+
+  MatrixXd R(F.rows(),6);
+  for (unsigned i=0; i<F.rows(); ++i)
+  {
+    RowVectorXd v;
+    canonical2frame_polar(TPs[i],thetas(i),S_polar.row(i),v);
+    R.row(i) = v;
+  }
+  return R;
+}
+
+  void FrameInterpolator::PolarDecomposition(Eigen::MatrixXd V, Eigen::MatrixXd& U, Eigen::MatrixXd& P)
+{
+  using namespace std;
+  using namespace Eigen;
+
+  // Polar Decomposition
+  JacobiSVD<MatrixXd> svd(V,Eigen::ComputeFullU | Eigen::ComputeFullV);
+
+  U = svd.matrixU() * svd.matrixV().transpose();
+  P = svd.matrixV() * svd.singularValues().asDiagonal() * svd.matrixV().transpose();
+
+  // If not SPD, change sign of both
+  // Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> S_es(P);
+  // if (S_es.eigenvalues().minCoeff() < 0)
+  // {
+  //   U = -U;
+  //   P = -P;
+  // }
+}
+
+}
+
+
+IGL_INLINE void igl::frame_field(
+                                 const Eigen::MatrixXd& V,
+                                 const Eigen::MatrixXi& F,
+                                 const Eigen::VectorXi& b,
+                                 const Eigen::MatrixXd& bc1,
+                                 const Eigen::MatrixXd& bc2,
+                                 Eigen::MatrixXd& F1,
+                                 Eigen::MatrixXd& F2,
+                                 Eigen::VectorXd& S
+                                 )
+
+{
+  using namespace std;
+  using namespace Eigen;
+
+  assert(b.size() > 0);
+  
+  // Init Solver
+  FrameInterpolator field(V,F);
+  
+  for (unsigned i=0; i<b.size(); ++i)
+  {
+    VectorXd t(6); t << bc1.row(i).transpose(), bc2.row(i).transpose();
+    field.setConstraint_polar(b(i), t,FrameInterpolator::SOFT);
+  }
+  
+  // Solve
+  field.solve_polar(false,true);
+  
+  // Copy back
+  MatrixXd R = field.getFieldPerFace_polar();
+  F1 = R.block(0, 0, R.rows(), 3);
+  F2 = R.block(0, 3, R.rows(), 3);
+}
+

include/igl/comiso/frame_field.h

@@ -0,0 +1,51 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2014 Daniele Panozzo <daniele.panozzo@gmail.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_FRAMEFIELD_H
+#define IGL_FRAMEFIELD_H
+
+#include <Eigen/Dense>
+#include <vector>
+
+namespace igl
+{
+// Generate a piecewise-constant frame-field field from a sparse set of constraints on faces
+// using the algorithm proposed in:
+// Frame Fields: Anisotropic and Non-Orthogonal Cross Fields
+// Daniele Panozzo, Enrico Puppo, Marco Tarini, Olga Sorkine-Hornung,
+// ACM Transactions on Graphics (SIGGRAPH, 2014)
+//
+// Inputs:
+//   V       #V by 3 list of mesh vertex coordinates
+//   F       #F by 3 list of mesh faces (must be triangles)
+//   b       #B by 1 list of constrained face indices
+//   bc1     #B by 3 list of the constrained first representative vector of the frame field (up to permutation and sign)
+//   bc2     #B by 3 list of the constrained second representative vector of the frame field (up to permutation and sign)
+//
+// Outputs:
+//   F1      #F by 3 the first representative vector of the frame field (up to permutation and sign)
+//   F2      #F by 3 the second representative vector of the frame field (up to permutation and sign)
+//   S       #V by 1 the singularity index for each vertex (0 = regular)
+//
+// TODO: it now supports only soft constraints, should be extended to support both hard and soft constraints
+IGL_INLINE void frame_field(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  const Eigen::VectorXi& b,
+  const Eigen::MatrixXd& bc1,
+  const Eigen::MatrixXd& bc2,
+  Eigen::MatrixXd& F1,
+  Eigen::MatrixXd& F2,
+  Eigen::VectorXd& S
+  );
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "frame_field.cpp"
+#endif
+
+#endif

include/igl/comiso/nrosy.h

@@ -30,7 +30,7 @@ namespace igl
 //           (0 -> smoothness only, 1->constraints only)
 
 // Outputs:
-//   R       #V by 3 the representative vectors of the interpolated field
+//   R       #F by 3 the representative vectors of the interpolated field
 //   S       #V by 1 the singularity index for each vertex (0 = regular)
 
 IGL_INLINE void nrosy(

include/igl/frame_field_deformer.cpp

@@ -0,0 +1,411 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2014 Daniele Panozzo <daniele.panozzo@gmail.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 "frame_field_deformer.h"
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+#include <vector>
+
+#include <igl/cotmatrix.h>
+#include <igl/vf.h>
+#include <igl/tt.h>
+
+namespace igl
+{
+
+class Frame_field_deformer
+{
+public:
+
+  IGL_INLINE Frame_field_deformer();
+  IGL_INLINE ~Frame_field_deformer();
+
+  // Initialize the optimizer
+  IGL_INLINE void init(const Eigen::MatrixXd& _V, const Eigen::MatrixXi& _F, const Eigen::MatrixXd& _D1, const Eigen::MatrixXd& _D2, double _Lambda, double _perturb_rotations, int _fixed = 1);
+
+  // Run N optimization steps
+  IGL_INLINE void optimize(int N, bool reset = false);
+
+  // Reset optimization
+  IGL_INLINE void reset_opt();
+
+  // Precomputation of all components
+  IGL_INLINE void precompute_opt();
+
+  // Precomputation for deformation energy
+  IGL_INLINE void precompute_ARAP(Eigen::SparseMatrix<double> & Lff, Eigen::MatrixXd & LfcVc);
+
+  // Precomputation for regularization
+  IGL_INLINE void precompute_SMOOTH(Eigen::SparseMatrix<double> & MS, Eigen::MatrixXd & bS);
+
+  // extracts a r x c block from sparse matrix mat into sparse matrix m1
+  // (r0,c0) is upper left entry of block
+  IGL_INLINE void extractBlock(Eigen::SparseMatrix<double> & mat, int r0, int c0, int r, int c, Eigen::SparseMatrix<double> & m1);
+
+  // computes optimal rotations for faces of m wrt current coords in mw.V
+  // returns a 3x3 matrix
+  IGL_INLINE void compute_optimal_rotations();
+
+  // global optimization step - linear system
+  IGL_INLINE void compute_optimal_positions();
+
+  // compute the output XField from deformation gradient
+  IGL_INLINE void computeXField(std::vector< Eigen::Matrix<double,3,2> > & XF);
+
+  // computes in WW the ideal warp at each tri to make the frame field a cross
+  IGL_INLINE void compute_idealWarp(std::vector< Eigen::Matrix<double,3,3> > & WW);
+
+  // -------------------------------- Variables ----------------------------------------------------
+
+  // Mesh I/O:
+
+  Eigen::MatrixXd V;                         // Original mesh - vertices
+  Eigen::MatrixXi F;                         // Original mesh - faces
+
+  std::vector<std::vector<int> > VT;                   // Vertex to triangle topology
+  std::vector<std::vector<int> > VTi;                  // Vertex to triangle topology
+
+  Eigen::MatrixXd V_w;                       // Warped mesh - vertices
+
+  std::vector< Eigen::Matrix<double,3,2> > FF;  	// frame field FF in 3D (parallel to m.F)
+  std::vector< Eigen::Matrix<double,3,3> > WW;    // warping matrices to make a cross field (parallel to m.F)
+  std::vector< Eigen::Matrix<double,3,2> > XF;  	// pseudo-cross field from solution (parallel to m.F)
+
+  int fixed;
+
+  double perturb_rotations; // perturbation to rotation matrices
+
+  // Numerics
+  int nfree,nconst;					          // number of free/constrained vertices in the mesh - default all-but-1/1
+  Eigen::MatrixXd C;							            // cotangent matrix of m
+  Eigen::SparseMatrix<double> L;					          // Laplacian matrix of m
+
+  Eigen::SparseMatrix<double> M;					          // matrix for global optimization - pre-conditioned
+  Eigen::MatrixXd RHS;						            // pre-computed part of known term in global optimization
+  std::vector< Eigen::Matrix<double,3,3> > RW;    // optimal rotation-warping matrices (parallel to m.F) -- INCORPORATES WW
+  Eigen::SimplicialCholesky<Eigen::SparseMatrix<double> > solver;   // solver for linear system in global opt.
+
+  // Parameters
+private:
+  double Lambda = 0.1;				        // weight of energy regularization
+
+};
+
+  IGL_INLINE Frame_field_deformer::Frame_field_deformer() {}
+
+  IGL_INLINE Frame_field_deformer::~Frame_field_deformer() {}
+
+  IGL_INLINE void Frame_field_deformer::init(const Eigen::MatrixXd& _V,
+                          const Eigen::MatrixXi& _F,
+                          const Eigen::MatrixXd& _D1,
+                          const Eigen::MatrixXd& _D2,
+                          double _Lambda,
+                          double _perturb_rotations,
+                          int _fixed)
+{
+  V = _V;
+  F = _F;
+
+  assert(_D1.rows() == _D2.rows());
+
+  FF.clear();
+  for (unsigned i=0; i < _D1.rows(); ++i)
+  {
+    Eigen::Matrix<double,3,2> ff;
+    ff.col(0) = _D1.row(i);
+    ff.col(1) = _D2.row(i);
+    FF.push_back(ff);
+  }
+
+  fixed = _fixed;
+  Lambda = _Lambda;
+  perturb_rotations = _perturb_rotations;
+
+  reset_opt();
+  precompute_opt();
+}
+
+
+IGL_INLINE void Frame_field_deformer::optimize(int N, bool reset)
+{
+  //Reset optimization
+	if (reset)
+    reset_opt();
+
+	// Iterative Local/Global optimization
+  for (int i=0; i<N;i++)
+  {
+    compute_optimal_rotations();
+    compute_optimal_positions();
+		computeXField(XF);
+  }
+}
+
+IGL_INLINE void Frame_field_deformer::reset_opt()
+{
+  V_w = V;
+
+  for (unsigned i=0; i<V_w.rows(); ++i)
+    for (unsigned j=0; j<V_w.cols(); ++j)
+      V_w(i,j) += (double(rand())/double(RAND_MAX))*10e-4*perturb_rotations;
+
+}
+
+// precomputation of all components
+IGL_INLINE void Frame_field_deformer::precompute_opt()
+{
+  using namespace Eigen;
+	nfree = V.rows() - fixed;						    // free vertices (at the beginning ov m.V) - global
+  nconst = V.rows()-nfree;						// #constrained vertices
+  igl::vf(V,F,VT,VTi);                // compute vertex to face relationship
+
+  igl::cotangent(V,F,C);							     // cotangent matrix for opt. rotations - global
+
+  igl::cotmatrix(V,F,L);
+
+	SparseMatrix<double> MA;						// internal matrix for ARAP-warping energy
+	MatrixXd LfcVc;										  // RHS (partial) for ARAP-warping energy
+	SparseMatrix<double> MS;						// internal matrix for smoothing energy
+	MatrixXd bS;										    // RHS (full) for smoothing energy
+
+	precompute_ARAP(MA,LfcVc);					// precompute terms for the ARAP-warp part
+	precompute_SMOOTH(MS,bS);						// precompute terms for the smoothing part
+	compute_idealWarp(WW);              // computes the ideal warps
+  RW.resize(F.rows());								// init rotation matrices - global
+
+  M =	  (1-Lambda)*MA + Lambda*MS;		// matrix for linear system - global
+
+	RHS = (1-Lambda)*LfcVc + Lambda*bS;	// RHS (partial) for linear system - global
+  solver.compute(M);									// system pre-conditioning
+  if (solver.info()!=Eigen::Success) {fprintf(stderr,"Decomposition failed in pre-conditioning!\n"); exit(-1);}
+
+	fprintf(stdout,"Preconditioning done.\n");
+
+}
+
+IGL_INLINE void Frame_field_deformer::precompute_ARAP(Eigen::SparseMatrix<double> & Lff, Eigen::MatrixXd & LfcVc)
+{
+  using namespace Eigen;
+	fprintf(stdout,"Precomputing ARAP terms\n");
+	SparseMatrix<double> LL = -4*L;
+	Lff = SparseMatrix<double>(nfree,nfree);
+  extractBlock(LL,0,0,nfree,nfree,Lff);
+	SparseMatrix<double> Lfc = SparseMatrix<double>(nfree,nconst);
+  extractBlock(LL,0,nfree,nfree,nconst,Lfc);
+	LfcVc = - Lfc * V_w.block(nfree,0,nconst,3);
+}
+
+IGL_INLINE void Frame_field_deformer::precompute_SMOOTH(Eigen::SparseMatrix<double> & MS, Eigen::MatrixXd & bS)
+{
+  using namespace Eigen;
+	fprintf(stdout,"Precomputing SMOOTH terms\n");
+
+	SparseMatrix<double> LL = 4*L*L;
+
+  // top-left
+	MS = SparseMatrix<double>(nfree,nfree);
+  extractBlock(LL,0,0,nfree,nfree,MS);
+
+  // top-right
+	SparseMatrix<double> Mfc = SparseMatrix<double>(nfree,nconst);
+  extractBlock(LL,0,nfree,nfree,nconst,Mfc);
+
+	MatrixXd MfcVc = Mfc * V_w.block(nfree,0,nconst,3);
+	bS = (LL*V).block(0,0,nfree,3)-MfcVc;
+
+}
+
+  IGL_INLINE void Frame_field_deformer::extractBlock(Eigen::SparseMatrix<double> & mat, int r0, int c0, int r, int c, Eigen::SparseMatrix<double> & m1)
+{
+  std::vector<Eigen::Triplet<double> > tripletList;
+  for (int k=c0; k<c0+c; ++k)
+    for (Eigen::SparseMatrix<double>::InnerIterator it(mat,k); it; ++it)
+    {
+      if (it.row()>=r0 && it.row()<r0+r)
+        tripletList.push_back(Eigen::Triplet<double>(it.row()-r0,it.col()-c0,it.value()));
+    }
+  m1.setFromTriplets(tripletList.begin(), tripletList.end());
+}
+
+IGL_INLINE void Frame_field_deformer::compute_optimal_rotations()
+{
+  using namespace Eigen;
+  Matrix<double,3,3> r,S,P,PP,D;
+
+  for (int i=0;i<F.rows();i++)
+	{
+		// input tri --- could be done once and saved in a matrix
+		P.col(0) = (V.row(F(i,1))-V.row(F(i,0))).transpose();
+		P.col(1) = (V.row(F(i,2))-V.row(F(i,1))).transpose();
+		P.col(2) = (V.row(F(i,0))-V.row(F(i,2))).transpose();
+
+		P = WW[i] * P;		// apply ideal warp
+
+		// current tri
+		PP.col(0) = (V_w.row(F(i,1))-V_w.row(F(i,0))).transpose();
+		PP.col(1) = (V_w.row(F(i,2))-V_w.row(F(i,1))).transpose();
+		PP.col(2) = (V_w.row(F(i,0))-V_w.row(F(i,2))).transpose();
+
+		// cotangents
+		D <<    C(i,2), 0,      0,
+    0,      C(i,0), 0,
+    0,      0,      C(i,1);
+
+		S = PP*D*P.transpose();
+		Eigen::JacobiSVD<Matrix<double,3,3> > svd(S, Eigen::ComputeFullU | Eigen::ComputeFullV );
+		Matrix<double,3,3>  su = svd.matrixU();
+		Matrix<double,3,3>  sv = svd.matrixV();
+		r = su*sv.transpose();
+
+		if (r.determinant()<0)  // correct reflections
+		{
+			su(0,2)=-su(0,2); su(1,2)=-su(1,2); su(2,2)=-su(2,2);
+			r = su*sv.transpose();
+		}
+		RW[i] = r*WW[i];		// RW INCORPORATES IDEAL WARP WW!!!
+	}
+}
+
+IGL_INLINE void Frame_field_deformer::compute_optimal_positions()
+{
+  using namespace Eigen;
+	// compute variable RHS of ARAP-warp part of the system
+  MatrixXd b(nfree,3);          // fx3 known term of the system
+	MatrixXd X;										// result
+  int t;		  									// triangles incident to edge (i,j)
+	int vi,i1,i2;									// index of vertex i wrt tri t0
+
+  for (int i=0;i<nfree;i++)
+  {
+    b.row(i) << 0.0, 0.0, 0.0;
+    for (int k=0;k<VT[i].size();k++)					// for all incident triangles
+    {
+      t = VT[i][k];												// incident tri
+			vi = (i==F(t,0))?0:(i==F(t,1))?1:(i==F(t,2))?2:3;	// index of i in t
+			assert(vi!=3);
+			i1 = F(t,(vi+1)%3);
+			i2 = F(t,(vi+2)%3);
+			b.row(i)+=(C(t,(vi+2)%3)*RW[t]*(V.row(i1)-V.row(i)).transpose()).transpose();
+			b.row(i)+=(C(t,(vi+1)%3)*RW[t]*(V.row(i2)-V.row(i)).transpose()).transpose();
+    }
+  }
+  b/=2.0;
+	b=-4*b;
+
+	b*=(1-Lambda);		// blend
+
+  b+=RHS;				// complete known term
+
+	X = solver.solve(b);
+	if (solver.info()!=Eigen::Success) {printf("Solving linear system failed!\n"); return;}
+
+	// copy result to mw.V
+  for (int i=0;i<nfree;i++)
+    V_w.row(i)=X.row(i);
+
+}
+
+  IGL_INLINE void Frame_field_deformer::computeXField(std::vector< Eigen::Matrix<double,3,2> > & XF)
+{
+  using namespace Eigen;
+  Matrix<double,3,3> P,PP,DG;
+	XF.resize(F.rows());
+
+  for (int i=0;i<F.rows();i++)
+	{
+		int i0,i1,i2;
+		// indexes of vertices of face i
+		i0 = F(i,0); i1 = F(i,1); i2 = F(i,2);
+
+		// input frame
+		P.col(0) = (V.row(i1)-V.row(i0)).transpose();
+		P.col(1) = (V.row(i2)-V.row(i0)).transpose();
+		P.col(2) = P.col(0).cross(P.col(1));
+    
+		// output triangle brought to origin
+		PP.col(0) = (V_w.row(i1)-V_w.row(i0)).transpose();
+		PP.col(1) = (V_w.row(i2)-V_w.row(i0)).transpose();
+		PP.col(2) = PP.col(0).cross(PP.col(1));
+
+		// deformation gradient
+		DG = PP * P.inverse();
+		XF[i] = DG * FF[i];
+	}
+}
+
+// computes in WW the ideal warp at each tri to make the frame field a cross
+  IGL_INLINE void Frame_field_deformer::compute_idealWarp(std::vector< Eigen::Matrix<double,3,3> > & WW)
+{
+  using namespace Eigen;
+  
+  WW.resize(F.rows());
+	for (int i=0;i<FF.size();i++)
+	{
+		Vector3d v0,v1,v2;
+		v0 = FF[i].col(0);
+		v1 = FF[i].col(1);
+		v2=v0.cross(v1); v2.normalize();			// normal
+
+		Matrix3d A,AI;								// compute affine map A that brings:
+		A <<    v0[0], v1[0], v2[0],				//	first vector of FF to x unary vector
+    v0[1], v1[1], v2[1],				//	second vector of FF to xy plane
+    v0[2], v1[2], v2[2];				//	triangle normal to z unary vector
+		AI = A.inverse();
+
+		// polar decomposition to discard rotational component (unnecessary but makes it easier)
+		Eigen::JacobiSVD<Matrix<double,3,3> > svd(AI, Eigen::ComputeFullU | Eigen::ComputeFullV );
+		Matrix<double,3,3>  au = svd.matrixU();
+		Matrix<double,3,3>  av = svd.matrixV();
+		DiagonalMatrix<double,3>	as(svd.singularValues());
+		WW[i] = av*as*av.transpose();
+	}
+}
+
+}
+
+
+IGL_INLINE void igl::frame_field_deformer(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  const Eigen::MatrixXd& FF1,
+  const Eigen::MatrixXd& FF2,
+  Eigen::MatrixXd&       V_d,
+  Eigen::MatrixXd&       FF1_d,
+  Eigen::MatrixXd&       FF2_d,
+  const int              iterations,
+  const double           lambda,
+  const bool             perturb_initial_guess)
+{
+  using namespace Eigen;
+  // Solvers
+  Frame_field_deformer deformer;
+
+  // Init optimizer
+  deformer.init(V, F, FF1, FF2, lambda, perturb_initial_guess ? 0.1 : 0);
+
+  // Optimize
+  deformer.optimize(iterations,true);
+
+  // Copy positions
+  V_d = deformer.V_w;
+
+  // Allocate
+  FF1_d.resize(F.rows(),3);
+  FF2_d.resize(F.rows(),3);
+
+  // Copy frame field
+  for(unsigned i=0; i<deformer.XF.size(); ++i)
+  {
+    FF1_d.row(i) = deformer.XF[i].col(0);
+    FF2_d.row(i) = deformer.XF[i].col(1);
+  }
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

include/igl/frame_field_deformer.h

@@ -0,0 +1,49 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2014 Daniele Panozzo <daniele.panozzo@gmail.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_FRAME_FIELD_DEFORMER_H
+#define IGL_FRAME_FIELD_DEFORMER_H
+#include "igl_inline.h"
+
+#include <Eigen/Dense>
+
+namespace igl
+{
+  // Deform a mesh to transform the given per-face frame field to be as close as possible
+  // to a cross field, in the least square sense.
+  //
+  // Inputs:
+  //   V       #V by 3 coordinates of the vertices
+  //   F       #F by 3 list of mesh faces (must be triangles)
+  //   FF1     #F by 3 first representative vector of the frame field
+  //   FF2     #F by 3 second representative vector of the frame field
+  //   lambda  laplacian regularization parameter 0=no regularization 1=full regularization
+  //
+  // Outputs:
+  //   V_d     #F by 3 deformed, first representative vector
+  //   V_d     #F by 3 deformed, first representative vector
+  //   V_d     #F by 3 deformed, first representative vector
+  //
+  IGL_INLINE void frame_field_deformer(
+    const Eigen::MatrixXd& V,
+    const Eigen::MatrixXi& F,
+    const Eigen::MatrixXd& FF1,
+    const Eigen::MatrixXd& FF2,
+    Eigen::MatrixXd&       V_d,
+    Eigen::MatrixXd&       FF1_d,
+    Eigen::MatrixXd&       FF2_d,
+    const int              iterations = 50,
+    const double           lambda = 0.1,
+    const bool             perturb_initial_guess = true);
+
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "frame_field_deformer.cpp"
+#endif
+
+#endif

tutorial/505_MIQ/main.cpp

@@ -67,11 +67,11 @@ bool key_down(igl::Viewer& viewer, unsigned char key, int modifier)
 {
   if (key <'1' || key >'7')
     return false;
-  
+
   viewer.clear_mesh();
   viewer.options.show_lines = false;
   viewer.options.show_texture = false;
-  
+
   if (key == '1')
   {
     // Cross field
@@ -79,7 +79,7 @@ bool key_down(igl::Viewer& viewer, unsigned char key, int modifier)
     viewer.add_edges (B, B + global_scale*X1 ,Eigen::RowVector3d(1,0,0));
     viewer.add_edges (B, B + global_scale*X2 ,Eigen::RowVector3d(0,0,1));
   }
-  
+
   if (key == '2')
   {
     // Bisector field
@@ -87,7 +87,7 @@ bool key_down(igl::Viewer& viewer, unsigned char key, int modifier)
     viewer.add_edges (B, B + global_scale*BIS1 ,Eigen::RowVector3d(1,0,0));
     viewer.add_edges (B, B + global_scale*BIS2 ,Eigen::RowVector3d(0,0,1));
   }
-  
+
   if (key == '3')
   {
     // Bisector field combed
@@ -95,7 +95,7 @@ bool key_down(igl::Viewer& viewer, unsigned char key, int modifier)
     viewer.add_edges (B, B + global_scale*BIS1_combed ,Eigen::RowVector3d(1,0,0));
     viewer.add_edges (B, B + global_scale*BIS2_combed ,Eigen::RowVector3d(0,0,1));
   }
-  
+
   if (key == '4')
   {
     // Singularities and cuts
@@ -105,7 +105,7 @@ bool key_down(igl::Viewer& viewer, unsigned char key, int modifier)
     int l_count = Seams.sum();
     Eigen::MatrixXd P1(l_count,3);
     Eigen::MatrixXd P2(l_count,3);
-    
+
     for (unsigned i=0; i<Seams.rows(); ++i)
     {
       for (unsigned j=0; j<Seams.cols(); ++j)
@@ -120,7 +120,7 @@ bool key_down(igl::Viewer& viewer, unsigned char key, int modifier)
     }
 
     viewer.add_edges(P1, P2, Eigen::RowVector3d(1, 0, 0));
-    
+
     // Plot the singularities as colored dots (red for negative, blue for positive)
     for (unsigned i=0; i<singularityIndex.size();++i)
     {
@@ -129,9 +129,9 @@ bool key_down(igl::Viewer& viewer, unsigned char key, int modifier)
       else if (singularityIndex(i) > 2)
         viewer.add_points(V.row(i),Eigen::RowVector3d(0,1,0));
     }
-    
+
   }
-  
+
   if (key == '5')
   {
     // Singularities and cuts, original field
@@ -139,12 +139,12 @@ bool key_down(igl::Viewer& viewer, unsigned char key, int modifier)
     viewer.set_mesh(V, F);
     viewer.add_edges (B, B + global_scale*X1_combed ,Eigen::RowVector3d(1,0,0));
     viewer.add_edges (B, B + global_scale*X2_combed ,Eigen::RowVector3d(0,0,1));
-    
+
     // Plot cuts
     int l_count = Seams.sum();
     Eigen::MatrixXd P1(l_count,3);
     Eigen::MatrixXd P2(l_count,3);
-    
+
     for (unsigned i=0; i<Seams.rows(); ++i)
     {
       for (unsigned j=0; j<Seams.cols(); ++j)
@@ -157,9 +157,9 @@ bool key_down(igl::Viewer& viewer, unsigned char key, int modifier)
         }
       }
     }
-    
+
     viewer.add_edges(P1, P2, Eigen::RowVector3d(1, 0, 0));
-    
+
     // Plot the singularities as colored dots (red for negative, blue for positive)
     for (unsigned i=0; i<singularityIndex.size();++i)
     {
@@ -169,7 +169,7 @@ bool key_down(igl::Viewer& viewer, unsigned char key, int modifier)
         viewer.add_points(V.row(i),Eigen::RowVector3d(0,1,0));
     }
   }
-  
+
   if (key == '6')
   {
     // Global parametrization UV
@@ -188,6 +188,11 @@ bool key_down(igl::Viewer& viewer, unsigned char key, int modifier)
 
   viewer.set_colors(Eigen::RowVector3d(1,1,1));
 
+  // Replace the standard texture with an integer shift invariant texture
+  Eigen::Matrix<char,Eigen::Dynamic,Eigen::Dynamic> texture_R, texture_G, texture_B;
+  line_texture(texture_R, texture_G, texture_B);
+  viewer.set_texture(texture_R, texture_B, texture_G);
+
   return false;
 }
 
@@ -200,7 +205,7 @@ int main(int argc, char *argv[])
 
   // Compute face barycenters
   igl::barycenter(V, F, B);
-  
+
   // Compute scale for visualizing fields
   global_scale =  .5*igl::avg_edge_length(V, F);
 
@@ -266,14 +271,9 @@ int main(int argc, char *argv[])
   // Plot the mesh
   igl::Viewer viewer;
 
-  // Replace the standard texture with an integer shift invariant texture
-  Eigen::Matrix<char,Eigen::Dynamic,Eigen::Dynamic> texture_R, texture_G, texture_B;
-  line_texture(texture_R, texture_G, texture_B);
-  viewer.set_texture(texture_R, texture_B, texture_G);
-  
   // Plot the original mesh with a texture parametrization
   key_down(viewer,'7',0);
-  
+
   // Launch the viewer
   viewer.callback_key_down = &key_down;
   viewer.launch();

tutorial/506_FrameField/CMakeLists.txt

@@ -0,0 +1,15 @@
+cmake_minimum_required(VERSION 2.6)
+project(506_FrameField)
+
+include("../CMakeLists.shared")
+
+find_package(LIBCOMISO REQUIRED)
+
+include_directories( ${LIBCOMISO_INCLUDE_DIR} )
+
+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/506_FrameField/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 <igl/comiso/nrosy.h>
+#include <igl/comiso/miq.h>
+#include <igl/comiso/frame_field.h>
+#include <igl/frame_field_deformer.h>
+#include <igl/jet.h>
+
+// Input mesh
+Eigen::MatrixXd V;
+Eigen::MatrixXi F;
+
+// Face barycenters
+Eigen::MatrixXd B;
+
+// Scale for visualizing the fields
+double global_scale;
+
+// Input frame field constraints
+Eigen::VectorXi b;
+Eigen::MatrixXd bc1;
+Eigen::MatrixXd bc2;
+
+// Interpolated frame field
+Eigen::MatrixXd FF1, FF2;
+
+// Singularities od the frame field
+Eigen::VectorXd S;
+
+// Deformed mesh
+Eigen::MatrixXd V_deformed;
+Eigen::MatrixXd B_deformed;
+
+// Frame field on deformed
+Eigen::MatrixXd FF1_deformed;
+Eigen::MatrixXd FF2_deformed;
+
+// Cross field on deformed
+Eigen::MatrixXd X1_deformed;
+Eigen::MatrixXd X2_deformed;
+
+// Quad mesh on deformed
+Eigen::MatrixXd V_quad_deformed;
+Eigen::MatrixXi F_quad_deformed;
+
+// Quad mesh
+Eigen::MatrixXd V_quad;
+Eigen::MatrixXi F_quad;
+
+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 == '1')
+  {
+    // Frame field constraints
+    viewer.set_mesh(V, F);
+
+    MatrixXd F1_t = MatrixXd::Zero(FF1.rows(),FF1.cols());
+    MatrixXd F2_t = MatrixXd::Zero(FF2.rows(),FF2.cols());
+    // 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;
+      F1_t.row(b(i)) = FF1.row(b(i));
+      F2_t.row(b(i)) = FF2.row(b(i));
+    }
+
+    viewer.set_colors(C);
+
+    MatrixXd C1,C2;
+    VectorXd K1 = F1_t.rowwise().norm();
+    VectorXd K2 = F2_t.rowwise().norm();
+    igl::jet(K1,true,C1);
+    igl::jet(K2,true,C2);
+
+    viewer.add_edges (B - global_scale*F1_t, B + global_scale*F1_t ,C1);
+    viewer.add_edges (B - global_scale*F2_t, B + global_scale*F2_t ,C2);
+  }
+
+  if (key == '2')
+  {
+    // Frame field
+    viewer.set_mesh(V, F);
+    MatrixXd C1,C2;
+    VectorXd K1 = FF1.rowwise().norm();
+    VectorXd K2 = FF2.rowwise().norm();
+    igl::jet(K1,true,C1);
+    igl::jet(K2,true,C2);
+
+    viewer.add_edges (B - global_scale*FF1, B + global_scale*FF1 ,C1);
+    viewer.add_edges (B - global_scale*FF2, B + global_scale*FF2 ,C2);
+
+    // 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.set_colors(C);
+
+  }
+
+  if (key == '3')
+  {
+    // Deformed with frame field
+    viewer.set_mesh(V_deformed, F);
+    viewer.add_edges (B_deformed, B_deformed + global_scale*FF1_deformed ,Eigen::RowVector3d(1,0,0));
+    viewer.add_edges (B_deformed, B_deformed + global_scale*FF2_deformed ,Eigen::RowVector3d(0,0,1));
+    viewer.set_colors(RowVector3d(1,1,1));
+  }
+
+  if (key == '4')
+  {
+    // Deformed with cross field
+    viewer.set_mesh(V_deformed, F);
+    viewer.add_edges (B, B + global_scale*X1_deformed ,Eigen::RowVector3d(1,0,0));
+    viewer.add_edges (B, B + global_scale*X2_deformed ,Eigen::RowVector3d(0,0,1));
+  }
+
+  if (key == '5')
+  {
+    // Deformed with quad mesh
+    viewer.set_mesh(V_quad_deformed, F_quad_deformed);
+  }
+
+  if (key == '6')
+  {
+    // Deformed with quad mesh
+    viewer.set_mesh(V_quad, F_quad);
+  }
+
+  return false;
+}
+
+int main(int argc, char *argv[])
+{
+  using namespace Eigen;
+
+  // Load a mesh in OBJ format
+  igl::readOBJ("../shared/cube.obj", V, F);
+
+  // Compute face barycenters
+  igl::barycenter(V, F, B);
+
+  // Compute scale for visualizing fields
+  global_scale =  .5*igl::avg_edge_length(V, F);
+
+  // Load constraints
+  MatrixXd temp;
+  igl::readDMAT("../shared/cube.dmat",temp);
+
+  b   = temp.block(0,0,temp.rows(),1).cast<int>();
+  bc1 = temp.block(0,1,temp.rows(),3);
+  bc2 = temp.block(0,4,temp.rows(),3);
+
+  // Interpolate the frame field
+  igl::frame_field(V, F, b, bc1, bc2, FF1, FF2, S);
+
+  // Deform the mesh to transform the frame field in a cross field
+  igl::frame_field_deformer(V,F,FF1,FF2,V_deformed,FF1_deformed,FF2_deformed);
+
+  // Compute face barycenters deformed mesh
+  igl::barycenter(V_deformed, F, B_deformed);
+
+  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();
+}