diff --git a/include/igl/comiso/frame_field.cpp b/include/igl/comiso/frame_field.cpp
index 0bed30b3c..90eaf5769 100644
--- a/include/igl/comiso/frame_field.cpp
+++ b/include/igl/comiso/frame_field.cpp
@@ -12,55 +12,33 @@ 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();
+  IGL_INLINE void setConstraint(const int fid, const Eigen::VectorXd& v);
 
-  // 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);
+  IGL_INLINE void interpolateSymmetric();
 
   // Generate the frame field
-  IGL_INLINE void solve_polar(const bool bilaplacian,bool hard_const);
+  IGL_INLINE void solve();
 
   // 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);
+  IGL_INLINE void frame2canonical(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 void canonical2frame(const Eigen::MatrixXd& TP, const double theta, const Eigen::VectorXd& S, Eigen::RowVectorXd& v);
 
-  IGL_INLINE Eigen::MatrixXd getFieldPerFace_polar();
+  IGL_INLINE Eigen::MatrixXd getFieldPerFace();
 
   IGL_INLINE void PolarDecomposition(Eigen::MatrixXd V, Eigen::MatrixXd& U, Eigen::MatrixXd& P);
 
   // Symmetric
-  Eigen::MatrixXd S_polar;
-  std::vector<ConstraintType> S_polar_c;
+  Eigen::MatrixXd S;
+  std::vector<bool> S_c;
 
   // -------------------------------------------------
 
@@ -82,20 +60,8 @@ private:
   // 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();
+  IGL_INLINE void interpolateCross();
 
   // Compute difference between reference frames
   IGL_INLINE void computek();
@@ -105,18 +71,7 @@ private:
 
   // 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;
+  std::vector<bool> thetas_c;
 
   // Edge Topology
   Eigen::MatrixXi EV, FE, EF;
@@ -133,9 +88,6 @@ private:
   // Normals per face
   Eigen::MatrixXd N;
 
-  // Singularity index
-  Eigen::VectorXd singularityIndex;
-
   // Reference frame per triangle
   std::vector<Eigen::MatrixXd> TPs;
 
@@ -189,25 +141,11 @@ FrameInterpolator::FrameInterpolator(const Eigen::MatrixXd& _V, const Eigen::Mat
   // 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());
+  S = MatrixXd::Zero(F.rows(),3);
 
   compute_edge_consistency();
-
 }
 
 FrameInterpolator::~FrameInterpolator()
@@ -215,80 +153,6 @@ 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)
@@ -321,10 +185,6 @@ double FrameInterpolator::vector2theta(const Eigen::MatrixXd& TP, const Eigen::R
 
   // 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;
 }
 
@@ -334,310 +194,35 @@ Eigen::RowVectorXd FrameInterpolator::theta2vector(const Eigen::MatrixXd& TP, co
   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)
+void FrameInterpolator::interpolateCross()
 {
   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);
-    }
-  }
+    if(thetas_c[i])
+      nrosy.setConstraintHard(i,theta2vector(TPs[i],thetas(i)));
 
   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);
+  for (unsigned i=0; i<F.rows(); ++i)
+    thetas(i) = vector2theta(TPs[i],R.row(i));
 }
 
 void FrameInterpolator::resetConstraints()
 {
   thetas_c.resize(F.rows());
-  betas_c.resize(F.rows());
-  ds_c.resize(F.rows());
-  S_polar_c.resize(F.rows());
+  S_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;
+    thetas_c[i]  = false;
+    S_c[i] = false;
   }
 
 }
@@ -789,7 +374,7 @@ void FrameInterpolator::computek()
 }
 
 
-  void FrameInterpolator::frame2canonical_polar(const Eigen::MatrixXd& TP, const Eigen::RowVectorXd& v, double& theta, Eigen::VectorXd& S_v)
+  void FrameInterpolator::frame2canonical(const Eigen::MatrixXd& TP, const Eigen::RowVectorXd& v, double& theta, Eigen::VectorXd& S_v)
 {
   using namespace std;
   using namespace Eigen;
@@ -813,68 +398,44 @@ void FrameInterpolator::computek()
   // 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));
+  theta = 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)
+  void FrameInterpolator::canonical2frame(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);
+  MatrixXd S_temp(2,2);
+  S_temp << 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;
+  MatrixXd M = S_temp * R;
 
   Vector2d vp0(M(0,0),M(1,0));
   Vector2d vp1(M(0,1),M(1,1));
@@ -887,21 +448,18 @@ void FrameInterpolator::computek()
   v << v0, v1;
 }
 
-void FrameInterpolator::solve_polar(const bool bilaplacian, bool hard_const)
+void FrameInterpolator::solve()
 {
-    interpolateCross(hard_const);
-    //interpolateSymmetric_polar(bilaplacian);
-  interpolateSymmetric_polar(false);
+  interpolateCross();
+  interpolateSymmetric();
 }
 
-void FrameInterpolator::interpolateSymmetric_polar(const bool bilaplacian)
+void FrameInterpolator::interpolateSymmetric()
 {
   using namespace std;
   using namespace Eigen;
 
-  //compute_edge_consistency();
-
-  // Generate uniform Laplacian matrix (see notes)
+  // Generate uniform Laplacian matrix
   typedef Eigen::Triplet<double> triplet;
   std::vector<triplet> triplets;
 
@@ -939,79 +497,26 @@ void FrameInterpolator::interpolateSymmetric_polar(const bool bilaplacian)
         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 ));
+        // 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 ));
+        // 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 ));
-        }
+        // 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 ));
       }
     }
   }
@@ -1025,12 +530,12 @@ void FrameInterpolator::interpolateSymmetric_polar(const bool bilaplacian)
   double w = 100000;
   for (unsigned fid=0; fid < F.rows(); ++fid)
   {
-    if (S_polar_c[fid] != FREE)
+    if (S_c[fid])
     {
       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);
+        b(3*fid + i) += w*S(fid,i);
       }
     }
   }
@@ -1040,68 +545,38 @@ void FrameInterpolator::interpolateSymmetric_polar(const bool bilaplacian)
 
   SparseMatrix<double> M;
 
-  if (!bilaplacian)
-    M = L + soft;
-  else
-    M = L * L + soft;
+  M = 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;
+    std::cerr << "LU 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);
+  S = 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);
+    S.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)
+void FrameInterpolator::setConstraint(const int fid, const Eigen::VectorXd& v)
 {
   using namespace std;
   using namespace Eigen;
@@ -1109,26 +584,20 @@ void FrameInterpolator::setConstraint_polar(const int fid, const Eigen::VectorXd
   double   t_;
   VectorXd S_;
 
-  frame2canonical_polar(TPs[fid],v,t_,S_);
+  frame2canonical(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;
+  canonical2frame(TPs[fid], t_, S_, v2);
 
   thetas(fid)   = t_;
-  thetas_c[fid] = type;
+  thetas_c[fid] = true;
 
-  S_polar.row(fid) = S_;
-  S_polar_c[fid]   = type;
+  S.row(fid) = S_;
+  S_c[fid]   = true;
 
 }
 
-Eigen::MatrixXd FrameInterpolator::getFieldPerFace_polar()
+Eigen::MatrixXd FrameInterpolator::getFieldPerFace()
 {
   using namespace std;
   using namespace Eigen;
@@ -1137,7 +606,7 @@ Eigen::MatrixXd FrameInterpolator::getFieldPerFace_polar()
   for (unsigned i=0; i<F.rows(); ++i)
   {
     RowVectorXd v;
-    canonical2frame_polar(TPs[i],thetas(i),S_polar.row(i),v);
+    canonical2frame(TPs[i],thetas(i),S.row(i),v);
     R.row(i) = v;
   }
   return R;
@@ -1153,14 +622,6 @@ Eigen::MatrixXd FrameInterpolator::getFieldPerFace_polar()
 
   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;
-  // }
 }
 
 }
@@ -1173,8 +634,7 @@ IGL_INLINE void igl::frame_field(
                                  const Eigen::MatrixXd& bc1,
                                  const Eigen::MatrixXd& bc2,
                                  Eigen::MatrixXd& F1,
-                                 Eigen::MatrixXd& F2,
-                                 Eigen::VectorXd& S
+                                 Eigen::MatrixXd& F2
                                  )
 
 {
@@ -1182,22 +642,21 @@ IGL_INLINE void igl::frame_field(
   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);
+    field.setConstraint(b(i), t);
   }
-  
+
   // Solve
-  field.solve_polar(false,true);
-  
+  field.solve();
+
   // Copy back
-  MatrixXd R = field.getFieldPerFace_polar();
+  MatrixXd R = field.getFieldPerFace();
   F1 = R.block(0, 0, R.rows(), 3);
   F2 = R.block(0, 3, R.rows(), 3);
 }
-
diff --git a/include/igl/comiso/frame_field.h b/include/igl/comiso/frame_field.h
index f409813c2..ae9470935 100644
--- a/include/igl/comiso/frame_field.h
+++ b/include/igl/comiso/frame_field.h
@@ -29,7 +29,6 @@ namespace igl
 // 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(
@@ -39,8 +38,7 @@ IGL_INLINE void frame_field(
   const Eigen::MatrixXd& bc1,
   const Eigen::MatrixXd& bc2,
   Eigen::MatrixXd& F1,
-  Eigen::MatrixXd& F2,
-  Eigen::VectorXd& S
+  Eigen::MatrixXd& F2
   );
 }
 
diff --git a/include/igl/rotate_vectors.cpp b/include/igl/rotate_vectors.cpp
index b193d8acf..8f8c09fb2 100644
--- a/include/igl/rotate_vectors.cpp
+++ b/include/igl/rotate_vectors.cpp
@@ -21,7 +21,7 @@ IGL_INLINE Eigen::MatrixXd igl::rotate_vectors(
     double a = atan2(B2.row(i).dot(V.row(i)),B1.row(i).dot(V.row(i)));
 
     // rotate
-    a += A.size() == 1 ? A(0) : A(1);
+    a += (A.size() == 1) ? A(0) : A(i);
 
     // move it back to global coordinates
     RV.row(i) = cos(a) * B1.row(i) + sin(a) * B2.row(i);
diff --git a/tutorial/506_FrameField/main.cpp b/tutorial/506_FrameField/main.cpp
index f718f4d2c..b7129f5ad 100755
--- a/tutorial/506_FrameField/main.cpp
+++ b/tutorial/506_FrameField/main.cpp
@@ -27,9 +27,6 @@ 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;
@@ -74,8 +71,8 @@ bool key_down(igl::Viewer& viewer, unsigned char key, int modifier)
     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));
+      F1_t.row(b(i)) = bc1.row(i);
+      F2_t.row(b(i)) = bc2.row(i);
     }
 
     viewer.set_colors(C);
@@ -154,7 +151,7 @@ int main(int argc, char *argv[])
   igl::barycenter(V, F, B);
 
   // Compute scale for visualizing fields
-  global_scale =  .5*igl::avg_edge_length(V, F);
+  global_scale =  .2*igl::avg_edge_length(V, F);
 
   // Load constraints
   MatrixXd temp;
@@ -164,8 +161,25 @@ int main(int argc, char *argv[])
   bc1 = temp.block(0,1,temp.rows(),3);
   bc2 = temp.block(0,4,temp.rows(),3);
 
+//  // Load a mesh in OBJ format
+//  igl::readOFF("../shared/planexy.off", V, F);
+//  
+//  // Compute face barycenters
+//  igl::barycenter(V, F, B);
+//  
+//  // Compute scale for visualizing fields
+//  global_scale =  .2*igl::avg_edge_length(V, F);
+//  
+//  b.resize(1);
+//  b << 0;
+//  bc1.resize(1,3);
+//  bc1 << (V.row(F(0,0)) - V.row(F(0,1))).normalized();
+//  MatrixXd B1,B2,B3;
+//  igl::local_basis(V,F,B1,B2,B3);
+//  bc2 = igl::rotate_vectors(bc1, VectorXd::Constant(1,M_PI/2), B1, B2);
+  
   // Interpolate the frame field
-  igl::frame_field(V, F, b, bc1, bc2, FF1, FF2, S);
+  igl::frame_field(V, F, b, bc1, bc2, FF1, FF2);
 
   // 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);