// This file is part of libigl, a simple c++ geometry processing library.
//
// Copyright (C) 2014 Daniele Panozzo <daniele.panozzo@gmail.com>, Olga Diamanti <olga.diam@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 <igl/comiso/miq.h>
#include <igl/local_basis.h>
#include <igl/triangle_triangle_adjacency.h>

// includes for VertexIndexing
#include <igl/HalfEdgeIterator.h>
#include <igl/is_border_vertex.h>
#include <igl/vertex_triangle_adjacency.h>


// includes for poissonSolver
#include <gmm/gmm.h>
#include <CoMISo/Solver/ConstrainedSolver.hh>
#include <CoMISo/Solver/MISolver.hh>
#include <CoMISo/Solver/GMM_Tools.hh>
#include <igl/doublearea.h>
#include <igl/per_face_normals.h>

//
#include <igl/cross_field_missmatch.h>
#include <igl/comb_frame_field.h>
#include <igl/cut_mesh_from_singularities.h>
#include <igl/find_cross_field_singularities.h>
#include <igl/compute_frame_field_bisectors.h>
#include <igl/rotate_vectors.h>

#define DEBUGPRINT 0


namespace igl {

  class SparseMatrixData{
  protected:
    unsigned int m_nrows;
    unsigned int m_ncols;
    std::vector<unsigned int> m_rowind;
    std::vector<unsigned int> m_colind;
    std::vector<double>       m_vals;

  public:
    unsigned int   nrows()    { return m_nrows      ; }
    unsigned int   ncols()    { return m_ncols      ; }
    unsigned int   nentries() { return m_vals.size(); }
    std::vector<unsigned int>&  rowind()   { return m_rowind     ; }
    std::vector<unsigned int>&  colind()   { return m_colind     ; }
    std::vector<double>&        vals()     { return m_vals       ; }

    // create an empty matrix with a fixed number of rows
    IGL_INLINE SparseMatrixData()
    {
      initialize(0,0);
    }

    // create an empty matrix with a fixed number of rows
    IGL_INLINE void initialize(int nr, int nc) {
      assert(nr >= 0 && nc >=0);
      m_nrows = nr;
      m_ncols = nc;

      m_rowind.resize(0);
      m_colind.resize(0);
      m_vals.resize(0);
    }

    // add a nonzero entry to the matrix
    // no checks are done for coinciding entries
    // the interpretation of the repeated entries (replace or add)
    // depends on how the actual sparse matrix datastructure is constructed

    IGL_INLINE void addEntryCmplx(unsigned int i, unsigned int j, std::complex<double> val) {
      m_rowind.push_back(2*i);   m_colind.push_back(2*j);   m_vals.push_back( val.real());
      m_rowind.push_back(2*i);   m_colind.push_back(2*j+1); m_vals.push_back(-val.imag());
      m_rowind.push_back(2*i+1); m_colind.push_back(2*j);   m_vals.push_back( val.imag());
      m_rowind.push_back(2*i+1); m_colind.push_back(2*j+1); m_vals.push_back( val.real());
    }

    IGL_INLINE void addEntryReal(unsigned int i, unsigned int j, double val) {
      m_rowind.push_back(i);   m_colind.push_back(j);   m_vals.push_back(val);
    }

    IGL_INLINE virtual ~SparseMatrixData() {
    }

  };

  // a small class to manage storage for matrix data
  // not using stl vectors: want to make all memory management
  // explicit to avoid hidden automatic reallocation
  // TODO: redo with STL vectors but with explicit mem. management

  class SparseSystemData {
  private:
    // matrix representation,  A[rowind[i],colind[i]] = vals[i]
    // right-hand side
    SparseMatrixData m_A;
    double       *m_b;
    double       *m_x;

  public:
    IGL_INLINE SparseMatrixData& A() { return m_A; }
    IGL_INLINE double*        b()        { return m_b       ; }
    IGL_INLINE double*        x()        { return m_x       ; }
    IGL_INLINE unsigned int   nrows()    { return  m_A.nrows(); }

  public:

    IGL_INLINE SparseSystemData(): m_A(), m_b(NULL), m_x(NULL){ }

    IGL_INLINE void initialize(unsigned int nr, unsigned int nc) {
      m_A.initialize(nr,nc);
      m_b      = new          double[nr];
      m_x      = new          double[nr];
      assert(m_b);
      std::fill( m_b,  m_b+nr, 0.);
    }

    IGL_INLINE void addRHSCmplx(unsigned int i, std::complex<double> val) {
      assert( 2*i+1 < m_A.nrows());
      m_b[2*i] += val.real(); m_b[2*i+1] += val.imag();
    }

    IGL_INLINE void setRHSCmplx(unsigned int i, std::complex<double> val) {
      assert( 2*i+1 < m_A.nrows());
      m_b[2*i] = val.real(); m_b[2*i+1] = val.imag();
    }

    IGL_INLINE std::complex<double> getRHSCmplx(unsigned int i) {
      assert( 2*i+1 < m_A.nrows());
      return std::complex<double>( m_b[2*i], m_b[2*i+1]);
    }

    IGL_INLINE double getRHSReal(unsigned int i) {
      assert( i < m_A.nrows());
      return m_b[i];
    }

    IGL_INLINE std::complex<double> getXCmplx(unsigned int i) {
      assert( 2*i+1 < m_A.nrows());
      return std::complex<double>( m_x[2*i], m_x[2*i+1]);
    }

    IGL_INLINE void cleanMem() {
      //m_A.cleanup();
      delete [] m_b;
      delete [] m_x;
    }

    IGL_INLINE virtual ~SparseSystemData() {
      delete [] m_b;
      delete [] m_x;
    }
  };

  struct SeamInfo
  {
    int v0,v0p,v1,v1p;
    int integerVar;
    unsigned char MMatch;

    IGL_INLINE SeamInfo(int _v0,
             int _v1,
             int _v0p,
             int _v1p,
             int _MMatch,
             int _integerVar);

    IGL_INLINE SeamInfo(const SeamInfo &S1);
  };

  struct MeshSystemInfo
  {
    ///total number of scalar variables
    int num_scalar_variables;
    ////number of vertices variables
    int num_vert_variables;
    ///num of integer for cuts
    int num_integer_cuts;
    ///this are used for drawing purposes
    std::vector<SeamInfo> EdgeSeamInfo;
#if 0
    ///this are values of integer variables after optimization
    std::vector<int> IntegerValues;
#endif
  };


  template <typename DerivedV, typename DerivedF>
  class VertexIndexing
  {
  public:
    // Input:
    const Eigen::PlainObjectBase<DerivedV> &V;
    const Eigen::PlainObjectBase<DerivedF> &F;
    const Eigen::PlainObjectBase<DerivedF> &TT;
    const Eigen::PlainObjectBase<DerivedF> &TTi;
    const Eigen::PlainObjectBase<DerivedV> &PD1;
    const Eigen::PlainObjectBase<DerivedV> &PD2;

    const Eigen::Matrix<int, Eigen::Dynamic, 3> &Handle_MMatch;
    const Eigen::Matrix<int, Eigen::Dynamic, 1> &Handle_Singular; // bool
    const Eigen::Matrix<int, Eigen::Dynamic, 1> &Handle_SingularDegree; // vertex;
    const Eigen::Matrix<int, Eigen::Dynamic, 3> &Handle_Seams; // 3 bool


    ///this handle for mesh TODO: move with the other global variables
    MeshSystemInfo Handle_SystemInfo;

    // Output:
    ///this maps the integer for edge - face
    Eigen::MatrixXi Handle_Integer; // TODO: remove it is useless

    ///per face indexes of vertex in the solver
    Eigen::MatrixXi HandleS_Index;

    ///per vertex variable indexes
    std::vector<std::vector<int> > HandleV_Integer;

    // internal
    std::vector<std::vector<int> > VF, VFi;
    std::vector<bool> V_border; // bool

    IGL_INLINE VertexIndexing(const Eigen::PlainObjectBase<DerivedV> &_V,
                   const Eigen::PlainObjectBase<DerivedF> &_F,
                   const Eigen::PlainObjectBase<DerivedF> &_TT,
                   const Eigen::PlainObjectBase<DerivedF> &_TTi,
                   const Eigen::PlainObjectBase<DerivedV> &_PD1,
                   const Eigen::PlainObjectBase<DerivedV> &_PD2,
                   const Eigen::Matrix<int, Eigen::Dynamic, 3> &_Handle_MMatch,
                   const Eigen::Matrix<int, Eigen::Dynamic, 1> &_Handle_Singular,
                   const Eigen::Matrix<int, Eigen::Dynamic, 1> &_Handle_SingularDegree,
                   const Eigen::Matrix<int, Eigen::Dynamic, 3> &_Handle_Seams
                   );

    ///vertex to variable mapping
    IGL_INLINE void InitMapping();

    IGL_INLINE void InitFaceIntegerVal();

    IGL_INLINE void InitSeamInfo();


  private:
    ///this maps back index to vertices
    std::vector<int> IndexToVert; // TODO remove it is useless

    ///this is used for drawing purposes
    std::vector<int> duplicated; // TODO remove it is useless

    IGL_INLINE void FirstPos(const int v, int &f, int &edge);

    IGL_INLINE int AddNewIndex(const int v0);

    IGL_INLINE bool HasIndex(int indexVert,int indexVar);

    IGL_INLINE void GetSeamInfo(const int f0,
                     const int f1,
                     const int indexE,
                     int &v0,int &v1,
                     int &v0p,int &v1p,
                     unsigned char &_MMatch,
                     int &integerVar);
    IGL_INLINE bool IsSeam(const int f0, const int f1);

    ///find initial position of the pos to
    // assing face to vert inxex correctly
    IGL_INLINE void FindInitialPos(const int vert, int &edge, int &face);


    ///intialize the mapping given an initial pos
    ///whih must be initialized with FindInitialPos
    IGL_INLINE void MapIndexes(const int  vert, const int edge_init, const int f_init);

    ///intialize the mapping for a given vertex
    IGL_INLINE void InitMappingSeam(const int vert);

    ///intialize the mapping for a given sampled mesh
    IGL_INLINE void InitMappingSeam();

    ///test consistency of face variables per vert mapping
    IGL_INLINE void TestSeamMappingFace(const int f);

    ///test consistency of face variables per vert mapping
    IGL_INLINE void TestSeamMappingVertex(int indexVert);

    ///check consistency of variable mapping across seams
    IGL_INLINE void TestSeamMapping();

  };


  template <typename DerivedV, typename DerivedF>
  class PoissonSolver
  {

  public:
    IGL_INLINE void SolvePoisson(Eigen::VectorXd Stiffness,
                      double vector_field_scale=0.1f,
                      double grid_res=1.f,
                      bool direct_round=true,
                      int localIter=0,
                      bool _integer_rounding=true,
                      std::vector<int> roundVertices = std::vector<int>(),
                      std::vector<std::vector<int> > hardFeatures = std::vector<std::vector<int> >());

    IGL_INLINE PoissonSolver(const Eigen::PlainObjectBase<DerivedV> &_V,
                  const Eigen::PlainObjectBase<DerivedF> &_F,
                  const Eigen::PlainObjectBase<DerivedF> &_TT,
                  const Eigen::PlainObjectBase<DerivedF> &_TTi,
                  const Eigen::PlainObjectBase<DerivedV> &_PD1,
                  const Eigen::PlainObjectBase<DerivedV> &_PD2,
                  const Eigen::MatrixXi &_HandleS_Index,
                  const Eigen::Matrix<int, Eigen::Dynamic, 1>&_Handle_Singular,
                  const MeshSystemInfo &_Handle_SystemInfo
                  );

    const Eigen::PlainObjectBase<DerivedV> &V;
    const Eigen::PlainObjectBase<DerivedF> &F;
    const Eigen::PlainObjectBase<DerivedF> &TT;
    const Eigen::PlainObjectBase<DerivedF> &TTi;
    const Eigen::PlainObjectBase<DerivedV> &PD1;
    const Eigen::PlainObjectBase<DerivedV> &PD2;
    const Eigen::Matrix<int, Eigen::Dynamic, 1> &Handle_Singular; // bool
    const Eigen::MatrixXi &HandleS_Index; //todo

    const MeshSystemInfo &Handle_SystemInfo;

    // Internal:
    Eigen::MatrixXd doublearea;
    Eigen::VectorXd Handle_Stiffness;
    Eigen::PlainObjectBase<DerivedV> N;
    std::vector<std::vector<int> > VF;
    std::vector<std::vector<int> > VFi;
    Eigen::MatrixXd UV; // this is probably useless

    // Output:
    // per wedge UV coordinates, 6 coordinates (1 face) per row
    Eigen::MatrixXd WUV;

    ///solver data
    SparseSystemData S;

    ///vector of unknowns
    std::vector< double > X;

    ////REAL PART
    ///number of fixed vertex
    unsigned int n_fixed_vars;

    ///the number of REAL variables for vertices
    unsigned int n_vert_vars;

    ///total number of variables of the system,
    ///do not consider constraints, but consider integer vars
    unsigned int num_total_vars;

    //////INTEGER PART
    ///the total number of integer variables
    unsigned int n_integer_vars;

    ///CONSTRAINT PART
    ///number of cuts constraints
    unsigned int num_cut_constraint;

    // number of user-defined constraints
    unsigned int num_userdefined_constraint;

    ///total number of constraints equations
    unsigned int num_constraint_equations;

    ///total size of the system including constraints
    unsigned int system_size;

    ///if you intend to make integer rotation
    ///and translations
    bool integer_jumps_bary;

    ///vector of blocked vertices
    std::vector<int> Hard_constraints;

    ///vector of indexes to round
    std::vector<int> ids_to_round;

    ///vector of indexes to round
    std::vector<std::vector<int > > userdefined_constraints;

    ///boolean that is true if rounding to integer is needed
    bool integer_rounding;

    ///START SYSTEM ACCESS METHODS
    ///add an entry to the LHS
    IGL_INLINE void AddValA(int Xindex,
                 int Yindex,
                 double val);

    ///add a complex entry to the LHS
    IGL_INLINE void AddComplexA(int VarXindex,
                     int VarYindex,
                     std::complex<double> val);

    ///add a velue to the RHS
    IGL_INLINE void AddValB(int Xindex,
                 double val);

    ///add the area term, scalefactor is used to sum up
    ///and normalize on the overlap zones
    IGL_INLINE void AddAreaTerm(int index[3][3][2],double ScaleFactor);

    ///set the diagonal of the matrix (which is zero at the beginning)
    ///such that the sum of a row or a colums is zero
    IGL_INLINE void SetDiagonal(double val[3][3]);

    ///given a vector of scalar values and
    ///a vector of indexes add such values
    ///as specified by the indexes
    IGL_INLINE void AddRHS(double b[6],
                int index[3]);

    ///add a 3x3 block matrix to the system matrix...
    ///indexes are specified in the 3x3 matrix of x,y pairs
    ///indexes must be multiplied by 2 cause u and v
    IGL_INLINE void Add33Block(double val[3][3], int index[3][3][2]);

    ///add a 3x3 block matrix to the system matrix...
    ///indexes are specified in the 3x3 matrix of x,y pairs
    ///indexes must be multiplied by 2 cause u and v
    IGL_INLINE void Add44Block(double val[4][4],int index[4][4][2]);
    ///END SYSTEM ACCESS METHODS

    ///START COMMON MATH FUNCTIONS
    ///return the complex encoding the rotation
    ///for a given missmatch interval
    IGL_INLINE std::complex<double> GetRotationComplex(int interval);
    ///END COMMON MATH FUNCTIONS


    ///START ENERGY MINIMIZATION PART
    ///initialize the LHS for a given face
    ///for minimization of Dirichlet's energy
    IGL_INLINE void perElementLHS(int f,
                       double val[3][3],
                       int index[3][3][2]);

    ///initialize the RHS for a given face
    ///for minimization of Dirichlet's energy
    IGL_INLINE void perElementRHS(int f,
                       double b[6],
                       double vector_field_scale=1);

    ///evaluate the LHS and RHS for a single face
    ///for minimization of Dirichlet's energy
    IGL_INLINE void PerElementSystemReal(int f,
                              double val[3][3],
                              int index[3][3][2],
                              double b[6],
                              double vector_field_scale=1.0);
    ///END ENERGY MINIMIZATION PART

    ///START FIXING VERTICES
    ///set a given vertex as fixed
    IGL_INLINE void AddFixedVertex(int v);

    ///find vertex to fix in case we're using
    ///a vector field NB: multiple components not handled
    IGL_INLINE void FindFixedVertField();

    ///find hard constraint depending if using or not
    ///a vector field
    IGL_INLINE void FindFixedVert();

    IGL_INLINE int GetFirstVertexIndex(int v);

    ///fix the vertices which are flagged as fixed
    IGL_INLINE void FixBlockedVertex();
    ///END FIXING VERTICES

    ///HANDLING SINGULARITY
    //set the singularity round to integer location
    IGL_INLINE void AddSingularityRound();

    IGL_INLINE void AddToRoundVertices(std::vector<int> ids);

    ///START GENERIC SYSTEM FUNCTIONS
    //build the laplacian matrix cyclyng over all rangemaps
    //and over all faces
    IGL_INLINE void BuildLaplacianMatrix(double vfscale=1);

    ///find different sized of the system
    IGL_INLINE void FindSizes();

    IGL_INLINE void AllocateSystem();

    ///intitialize the whole matrix
    IGL_INLINE void InitMatrix();

    ///map back coordinates after that
    ///the system has been solved
    IGL_INLINE void MapCoords();
    ///END GENERIC SYSTEM FUNCTIONS

    ///set the constraints for the inter-range cuts
    IGL_INLINE void BuildSeamConstraintsExplicitTranslation();

    ///set the constraints for the inter-range cuts
    IGL_INLINE void BuildUserDefinedConstraints();

    ///call of the mixed integer solver
    IGL_INLINE void MixedIntegerSolve(double cone_grid_res=1,
                           bool direct_round=true,
                           int localIter=0);

    IGL_INLINE void clearUserConstraint();

    IGL_INLINE void addSharpEdgeConstraint(int fid, int vid);

  };

  template <typename DerivedV, typename DerivedF, typename DerivedU>
  class MIQ_class
  {
  private:
    const Eigen::PlainObjectBase<DerivedV> &V;
    const Eigen::PlainObjectBase<DerivedF> &F;
    Eigen::MatrixXd WUV;
    // internal
    Eigen::PlainObjectBase<DerivedF> TT;
    Eigen::PlainObjectBase<DerivedF> TTi;

    // Stiffness per face
    Eigen::VectorXd Handle_Stiffness;
    Eigen::PlainObjectBase<DerivedV> B1, B2, B3;

  public:
    IGL_INLINE MIQ_class(const Eigen::PlainObjectBase<DerivedV> &V_,
              const Eigen::PlainObjectBase<DerivedF> &F_,
              const Eigen::PlainObjectBase<DerivedV> &PD1_combed,
              const Eigen::PlainObjectBase<DerivedV> &PD2_combed,
              const Eigen::PlainObjectBase<DerivedV> &BIS1_combed,
              const Eigen::PlainObjectBase<DerivedV> &BIS2_combed,
              const Eigen::Matrix<int, Eigen::Dynamic, 3> &Handle_MMatch,
              const Eigen::Matrix<int, Eigen::Dynamic, 1> &Handle_Singular,
              const Eigen::Matrix<int, Eigen::Dynamic, 1> &Handle_SingularDegree,
              const Eigen::Matrix<int, Eigen::Dynamic, 3> &Handle_Seams,
              Eigen::PlainObjectBase<DerivedU> &UV,
              Eigen::PlainObjectBase<DerivedF> &FUV,
              double GradientSize = 30.0,
              double Stiffness = 5.0,
              bool DirectRound = false,
              int iter = 5,
              int localIter = 5,
              bool DoRound = true,
              std::vector<int> roundVertices = std::vector<int>(),
              std::vector<std::vector<int> > hardFeatures = std::vector<std::vector<int> >());


    IGL_INLINE void extractUV(Eigen::PlainObjectBase<DerivedU> &UV_out,
                   Eigen::PlainObjectBase<DerivedF> &FUV_out);

  private:
    IGL_INLINE int NumFlips(const Eigen::MatrixXd& WUV);

    IGL_INLINE double Distortion(int f, double h, const Eigen::MatrixXd& WUV);

    IGL_INLINE double LaplaceDistortion(const int f, double h, const Eigen::MatrixXd& WUV);

    IGL_INLINE bool updateStiffeningJacobianDistorsion(double grad_size, const Eigen::MatrixXd& WUV);

    IGL_INLINE bool IsFlipped(const Eigen::Vector2d &uv0,
                          const Eigen::Vector2d &uv1,
                          const Eigen::Vector2d &uv2);

    IGL_INLINE bool IsFlipped(const int i, const Eigen::MatrixXd& WUV);

  };
};

IGL_INLINE igl::SeamInfo::SeamInfo(int _v0,
                        int _v1,
                        int _v0p,
                        int _v1p,
                        int _MMatch,
                        int _integerVar)
{
  v0=_v0;
  v1=_v1;
  v0p=_v0p;
  v1p=_v1p;
  integerVar=_integerVar;
  MMatch=_MMatch;
}

IGL_INLINE igl::SeamInfo::SeamInfo(const SeamInfo &S1)
{
  v0=S1.v0;
  v1=S1.v1;
  v0p=S1.v0p;
  v1p=S1.v1p;
  integerVar=S1.integerVar;
  MMatch=S1.MMatch;
}


template <typename DerivedV, typename DerivedF>
IGL_INLINE igl::VertexIndexing<DerivedV, DerivedF>::VertexIndexing(const Eigen::PlainObjectBase<DerivedV> &_V,
                                                        const Eigen::PlainObjectBase<DerivedF> &_F,
                                                        const Eigen::PlainObjectBase<DerivedF> &_TT,
                                                        const Eigen::PlainObjectBase<DerivedF> &_TTi,
                                                        const Eigen::PlainObjectBase<DerivedV> &_PD1,
                                                        const Eigen::PlainObjectBase<DerivedV> &_PD2,
                                                        const Eigen::Matrix<int, Eigen::Dynamic, 3> &_Handle_MMatch,
                                                        const Eigen::Matrix<int, Eigen::Dynamic, 1> &_Handle_Singular,
                                                        const Eigen::Matrix<int, Eigen::Dynamic, 1> &_Handle_SingularDegree,
                                                        const Eigen::Matrix<int, Eigen::Dynamic, 3> &_Handle_Seams

                                                        ):
V(_V),
F(_F),
TT(_TT),
TTi(_TTi),
PD1(_PD1),
PD2(_PD2),
Handle_MMatch(_Handle_MMatch),
Handle_Singular(_Handle_Singular),
Handle_SingularDegree(_Handle_SingularDegree),
Handle_Seams(_Handle_Seams)
{

  V_border = igl::is_border_vertex(V,F);
  igl::vertex_triangle_adjacency(V,F,VF,VFi);

  IndexToVert.clear();

  Handle_SystemInfo.num_scalar_variables=0;
  Handle_SystemInfo.num_vert_variables=0;
  Handle_SystemInfo.num_integer_cuts=0;

  duplicated.clear();

  HandleS_Index = Eigen::MatrixXi::Constant(F.rows(),3,-1);

  Handle_Integer = Eigen::MatrixXi::Constant(F.rows(),3,-1);

  HandleV_Integer.resize(V.rows());
}

template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::VertexIndexing<DerivedV, DerivedF>::FirstPos(const int v, int &f, int &edge)
{
  f    = VF[v][0];  // f=v->cVFp();
  edge = VFi[v][0]; // edge=v->cVFi();
}

template <typename DerivedV, typename DerivedF>
IGL_INLINE int igl::VertexIndexing<DerivedV, DerivedF>::AddNewIndex(const int v0)
{
  Handle_SystemInfo.num_scalar_variables++;
  HandleV_Integer[v0].push_back(Handle_SystemInfo.num_scalar_variables);
  IndexToVert.push_back(v0);
  return Handle_SystemInfo.num_scalar_variables;
}

template <typename DerivedV, typename DerivedF>
IGL_INLINE bool igl::VertexIndexing<DerivedV, DerivedF>::HasIndex(int indexVert,int indexVar)
{
  for (unsigned int i=0;i<HandleV_Integer[indexVert].size();i++)
    if (HandleV_Integer[indexVert][i]==indexVar)return true;
  return false;
}

template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::VertexIndexing<DerivedV, DerivedF>::GetSeamInfo(const int f0,
                                                          const int f1,
                                                          const int indexE,
                                                          int &v0,int &v1,
                                                          int &v0p,int &v1p,
                                                          unsigned char &_MMatch,
                                                          int &integerVar)
{
  int edgef0 = indexE;
  v0 = HandleS_Index(f0,edgef0);
  v1 = HandleS_Index(f0,(edgef0+1)%3);
  ////get the index on opposite side
  assert(TT(f0,edgef0) == f1);
  int edgef1 = TTi(f0,edgef0);
  v1p = HandleS_Index(f1,edgef1);
  v0p = HandleS_Index(f1,(edgef1+1)%3);

  integerVar = Handle_Integer(f0,edgef0);
  _MMatch = Handle_MMatch(f0,edgef0);
  assert(F(f0,edgef0)         == F(f1,((edgef1+1)%3)));
  assert(F(f0,((edgef0+1)%3)) == F(f1,edgef1));
}

template <typename DerivedV, typename DerivedF>
IGL_INLINE bool igl::VertexIndexing<DerivedV, DerivedF>::IsSeam(const int f0, const int f1)
{
  for (int i=0;i<3;i++)
  {
    int f_clos = TT(f0,i);

    if (f_clos == -1)
      continue; ///border

    if (f_clos == f1)
      return(Handle_Seams(f0,i));
  }
  assert(0);
  return false;
}

///find initial position of the pos to
// assing face to vert inxex correctly
template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::VertexIndexing<DerivedV, DerivedF>::FindInitialPos(const int vert,
                                                             int &edge,
                                                             int &face)
{
  int f_init;
  int edge_init;
  FirstPos(vert,f_init,edge_init); // todo manually IGL_INLINE the function
  igl::HalfEdgeIterator<DerivedF> VFI(&F,&TT,&TTi,f_init,edge_init);

  bool vertexB = V_border[vert];
  bool possible_split=false;
  bool complete_turn=false;
  do
  {
    int curr_f = VFI.Fi();
    int curr_edge=VFI.Ei();
    VFI.NextFE();
    int next_f=VFI.Fi();
    ///test if I've just crossed a border
    bool on_border=(TT(curr_f,curr_edge)==-1);
    //bool mismatch=false;
    bool seam=false;

    ///or if I've just crossed a seam
    ///if I'm on a border I MUST start from the one next t othe border
    if (!vertexB)
      //seam=curr_f->IsSeam(next_f);
      seam=IsSeam(curr_f,next_f);
    if (vertexB)
      assert(!Handle_Singular(vert));
    ;
    //assert(!vert->IsSingular());
    possible_split=((on_border)||(seam));
    complete_turn = next_f == f_init;
  } while ((!possible_split)&&(!complete_turn));
  face=VFI.Fi();
  edge=VFI.Ei();
  ///test that is not on a border
  //assert(face->FFp(edge)!=face);
}



///intialize the mapping given an initial pos
///whih must be initialized with FindInitialPos
template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::VertexIndexing<DerivedV, DerivedF>::MapIndexes(const int  vert,
                                                         const int edge_init,
                                                         const int f_init)
{
  ///check that is not on border..
  ///in such case maybe it's non manyfold
  ///insert an initial index
  int curr_index=AddNewIndex(vert);
  ///and initialize the jumping pos
  igl::HalfEdgeIterator<DerivedF> VFI(&F,&TT,&TTi,f_init,edge_init);
  bool complete_turn=false;
  do
  {
    int curr_f = VFI.Fi();
    int curr_edge = VFI.Ei();
    ///assing the current index
    HandleS_Index(curr_f,curr_edge) = curr_index;
    VFI.NextFE();
    int next_f = VFI.Fi();
    ///test if I've finiseh with the face exploration
    complete_turn = (next_f==f_init);
    ///or if I've just crossed a mismatch
    if (!complete_turn)
    {
      bool seam=false;
      //seam=curr_f->IsSeam(next_f);
      seam=IsSeam(curr_f,next_f);
      if (seam)
      {
        ///then add a new index
        curr_index=AddNewIndex(vert);
      }
    }
  } while (!complete_turn);
}

///intialize the mapping for a given vertex
template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::VertexIndexing<DerivedV, DerivedF>::InitMappingSeam(const int vert)
{
  ///first rotate until find the first pos after a mismatch
  ///or a border or return to the first position...
  int f_init = VF[vert][0];
  int indexE = VFi[vert][0];

  igl::HalfEdgeIterator<DerivedF> VFI(&F,&TT,&TTi,f_init,indexE);

  int edge_init;
  int face_init;
  FindInitialPos(vert,edge_init,face_init);
  MapIndexes(vert,edge_init,face_init);
}

///intialize the mapping for a given sampled mesh
template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::VertexIndexing<DerivedV, DerivedF>::InitMappingSeam()
{
  //num_scalar_variables=-1;
  Handle_SystemInfo.num_scalar_variables=-1;
  for (unsigned int i=0;i<V.rows();i++)
    InitMappingSeam(i);

  for (unsigned int j=0;j<V.rows();j++)
  {
    assert(HandleV_Integer[j].size()>0);
    if (HandleV_Integer[j].size()>1)
      duplicated.push_back(j);
  }
}

///test consistency of face variables per vert mapping
template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::VertexIndexing<DerivedV, DerivedF>::TestSeamMappingFace(const int f)
{
  for (int k=0;k<3;k++)
  {
    int indexV=HandleS_Index(f,k);
    int v = F(f,k);
    bool has_index=HasIndex(v,indexV);
    assert(has_index);
  }
}

///test consistency of face variables per vert mapping
template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::VertexIndexing<DerivedV, DerivedF>::TestSeamMappingVertex(int indexVert)
{
  for (unsigned int k=0;k<HandleV_Integer[indexVert].size();k++)
  {
    int indexV=HandleV_Integer[indexVert][k];

    ///get faces sharing vertex
    std::vector<int> faces = VF[indexVert];
    std::vector<int> indexes = VFi[indexVert];

    for (unsigned int j=0;j<faces.size();j++)
    {
      int f = faces[j];
      int index = indexes[j];
      assert(F(f,index) == indexVert);
      assert((index>=0)&&(index<3));

      if (HandleS_Index(f,index) == indexV)
        return;
    }
  }
  assert(0);
}


///check consistency of variable mapping across seams
template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::VertexIndexing<DerivedV, DerivedF>::TestSeamMapping()
{
  printf("\n TESTING SEAM INDEXES \n");
  ///test F-V mapping
  for (unsigned int j=0;j<F.rows();j++)
    TestSeamMappingFace(j);

  ///TEST  V-F MAPPING
  for (unsigned int j=0;j<V.rows();j++)
    TestSeamMappingVertex(j);

}


///vertex to variable mapping
template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::VertexIndexing<DerivedV, DerivedF>::InitMapping()
{
  //use_direction_field=_use_direction_field;

  IndexToVert.clear();
  duplicated.clear();

  InitMappingSeam();

  Handle_SystemInfo.num_vert_variables=Handle_SystemInfo.num_scalar_variables+1;

  ///end testing...
  TestSeamMapping();
}

template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::VertexIndexing<DerivedV, DerivedF>::InitFaceIntegerVal()
{
  Handle_SystemInfo.num_integer_cuts=0;
  for (unsigned int j=0;j<F.rows();j++)
  {
    for (int k=0;k<3;k++)
    {
      if (Handle_Seams(j,k))
      {
        Handle_Integer(j,k) = Handle_SystemInfo.num_integer_cuts;
        Handle_SystemInfo.num_integer_cuts++;
      }
      else
        Handle_Integer(j,k)=-1;
    }
  }
}


template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::VertexIndexing<DerivedV, DerivedF>::InitSeamInfo()
{
  Handle_SystemInfo.EdgeSeamInfo.clear();
  for (unsigned int f0=0;f0<F.rows();f0++)
  {
    for (int k=0;k<3;k++)
    {
      int f1 = TT(f0,k);

      if (f1 == -1)
        continue;

      bool seam = Handle_Seams(f0,k);
      if (seam)
      {
        int v0,v0p,v1,v1p;
        unsigned char MM;
        int integerVar;
        GetSeamInfo(f0,f1,k,v0,v1,v0p,v1p,MM,integerVar);
        Handle_SystemInfo.EdgeSeamInfo.push_back(SeamInfo(v0,v1,v0p,v1p,MM,integerVar));
      }
    }
  }
}



template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::SolvePoisson(Eigen::VectorXd Stiffness,
                                                          double vector_field_scale,
                                                          double grid_res,
                                                          bool direct_round,
                                                          int localIter,
                                                          bool _integer_rounding,
                                                          std::vector<int> roundVertices,
                                                          std::vector<std::vector<int> > hardFeatures)
{
  Handle_Stiffness = Stiffness;

  //initialization of flags and data structures
  integer_rounding=_integer_rounding;

  ids_to_round.clear();

  clearUserConstraint();
  // copy the user constraints number
  for (size_t i = 0; i < hardFeatures.size(); ++i)
  {
    addSharpEdgeConstraint(hardFeatures[i][0],hardFeatures[i][1]);
  }

  ///Initializing Matrix

  int t0=clock();

  ///initialize the matrix ALLOCATING SPACE
  InitMatrix();
  if (DEBUGPRINT)
    printf("\n ALLOCATED THE MATRIX \n");

  ///build the laplacian system
  BuildLaplacianMatrix(vector_field_scale);

  // add seam constraints
  BuildSeamConstraintsExplicitTranslation();

  // add user defined constraints
  BuildUserDefinedConstraints();

  ////add the lagrange multiplier
  FixBlockedVertex();

  if (DEBUGPRINT)
    printf("\n BUILT THE MATRIX \n");

  if (integer_rounding)
  {
    AddSingularityRound();
    AddToRoundVertices(roundVertices);
  }

  int t1=clock();
  if (DEBUGPRINT) printf("\n time:%d \n",t1-t0);
  if (DEBUGPRINT) printf("\n SOLVING \n");

  MixedIntegerSolve(grid_res,direct_round,localIter);

  int t2=clock();
  if (DEBUGPRINT) printf("\n time:%d \n",t2-t1);
  if (DEBUGPRINT) printf("\n ASSIGNING COORDS \n");

  MapCoords();

  int t3=clock();
  if (DEBUGPRINT) printf("\n time:%d \n",t3-t2);
  if (DEBUGPRINT) printf("\n FINISHED \n");
}

template <typename DerivedV, typename DerivedF>
IGL_INLINE igl::PoissonSolver<DerivedV, DerivedF>
::PoissonSolver(const Eigen::PlainObjectBase<DerivedV> &_V,
                const Eigen::PlainObjectBase<DerivedF> &_F,
                const Eigen::PlainObjectBase<DerivedF> &_TT,
                const Eigen::PlainObjectBase<DerivedF> &_TTi,
                const Eigen::PlainObjectBase<DerivedV> &_PD1,
                const Eigen::PlainObjectBase<DerivedV> &_PD2,
                const Eigen::MatrixXi &_HandleS_Index,
                const Eigen::Matrix<int, Eigen::Dynamic, 1>&_Handle_Singular,
                const MeshSystemInfo &_Handle_SystemInfo //todo: const?
):
V(_V),
F(_F),
TT(_TT),
TTi(_TTi),
PD1(_PD1),
PD2(_PD2),
HandleS_Index(_HandleS_Index),
Handle_Singular(_Handle_Singular),
Handle_SystemInfo(_Handle_SystemInfo)
{
  UV        = Eigen::MatrixXd(V.rows(),2);
  WUV       = Eigen::MatrixXd(F.rows(),6);
  igl::doublearea(V,F,doublearea);
  igl::per_face_normals(V,F,N);
  igl::vertex_triangle_adjacency(V,F,VF,VFi);
}


///START SYSTEM ACCESS METHODS
///add an entry to the LHS
template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::AddValA(int Xindex,
                                                     int Yindex,
                                                     double val)
{
  int size=(int)S.nrows();
  assert(0 <= Xindex && Xindex < size);
  assert(0 <= Yindex && Yindex < size);
  S.A().addEntryReal(Xindex,Yindex,val);
}

///add a complex entry to the LHS
template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::AddComplexA(int VarXindex,
                                                         int VarYindex,
                                                         std::complex<double> val)
{
  int size=(int)S.nrows()/2;
  assert(0 <= VarXindex && VarXindex < size);
  assert(0 <= VarYindex && VarYindex < size);
  S.A().addEntryCmplx(VarXindex,VarYindex,val);
}

///add a velue to the RHS
template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::AddValB(int Xindex,
                                                     double val)
{
  int size=(int)S.nrows();
  assert(0 <= Xindex && Xindex < size);
  S.b()[Xindex] += val;
}

///add the area term, scalefactor is used to sum up
///and normalize on the overlap zones
template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::AddAreaTerm(int index[3][3][2],double ScaleFactor)
{
  const double entry = 0.5*ScaleFactor;
  double val[3][3]= {
    {0,       entry, -entry},
    {-entry,      0,  entry},
    {entry,  -entry,      0}
  };

  for (int i=0;i<3;i++)
    for (int j=0;j<3;j++)
    {
      ///add for both u and v
      int Xindex=index[i][j][0]*2;
      int Yindex=index[i][j][1]*2;

      AddValA(Xindex+1,Yindex,-val[i][j]);
      AddValA(Xindex,Yindex+1,val[i][j]);
    }
}

///set the diagonal of the matrix (which is zero at the beginning)
///such that the sum of a row or a colums is zero
template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::SetDiagonal(double val[3][3])
{
  for (int i=0;i<3;i++)
  {
    double sum=0;
    for (int j=0;j<3;j++)
      sum+=val[i][j];
    val[i][i]=-sum;
  }
}

///given a vector of scalar values and
///a vector of indexes add such values
///as specified by the indexes
template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::AddRHS(double b[6],
                                                    int index[3])
{
  for (int i=0;i<3;i++)
  {
    double valU=b[i*2];
    double valV=b[(i*2)+1];
    AddValB((index[i]*2),valU);
    AddValB((index[i]*2)+1,valV);
  }
}

///add a 3x3 block matrix to the system matrix...
///indexes are specified in the 3x3 matrix of x,y pairs
///indexes must be multiplied by 2 cause u and v
template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::Add33Block(double val[3][3], int index[3][3][2])
{
  for (int i=0;i<3;i++)
    for (int j=0;j<3;j++)
    {
      ///add for both u and v
      int Xindex=index[i][j][0]*2;
      int Yindex=index[i][j][1]*2;
      assert((unsigned)Xindex<(n_vert_vars*2));
      assert((unsigned)Yindex<(n_vert_vars*2));
      AddValA(Xindex,Yindex,val[i][j]);
      AddValA(Xindex+1,Yindex+1,val[i][j]);
    }
}

///add a 3x3 block matrix to the system matrix...
///indexes are specified in the 3x3 matrix of x,y pairs
///indexes must be multiplied by 2 cause u and v
template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::Add44Block(double val[4][4],int index[4][4][2])
{
  for (int i=0;i<4;i++)
    for (int j=0;j<4;j++)
    {
      ///add for both u and v
      int Xindex=index[i][j][0]*2;
      int Yindex=index[i][j][1]*2;
      assert((unsigned)Xindex<(n_vert_vars*2));
      assert((unsigned)Yindex<(n_vert_vars*2));
      AddValA(Xindex,Yindex,val[i][j]);
      AddValA(Xindex+1,Yindex+1,val[i][j]);
    }
}
///END SYSTEM ACCESS METHODS

///START COMMON MATH FUNCTIONS
///return the complex encoding the rotation
///for a given missmatch interval
template <typename DerivedV, typename DerivedF>
IGL_INLINE std::complex<double> igl::PoissonSolver<DerivedV, DerivedF>::GetRotationComplex(int interval)
{
  assert((interval>=0)&&(interval<4));

  switch(interval)
  {
    case 0:return std::complex<double>(1,0);
    case 1:return std::complex<double>(0,1);
    case 2:return std::complex<double>(-1,0);
    default:return std::complex<double>(0,-1);
  }
}

///END COMMON MATH FUNCTIONS


///START ENERGY MINIMIZATION PART
///initialize the LHS for a given face
///for minimization of Dirichlet's energy
template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::perElementLHS(int f,
                                                           double val[3][3],
                                                           int index[3][3][2])
{
  ///initialize to zero
  for (int x=0;x<3;x++)
    for (int y=0;y<3;y++)
      val[x][y]=0;

  ///get the vertices
  int v[3];
  v[0] = F(f,0);
  v[1] = F(f,1);
  v[2] = F(f,2);

  ///get the indexes of vertex instance (to consider cuts)
  ///for the current face
  int Vindexes[3];
  Vindexes[0]=HandleS_Index(f,0);
  Vindexes[1]=HandleS_Index(f,1);
  Vindexes[2]=HandleS_Index(f,2);

  ///initialize the indexes for the block
  for (int x=0;x<3;x++)
    for (int y=0;y<3;y++)
    {
      index[x][y][0]=Vindexes[x];
      index[x][y][1]=Vindexes[y];
    }

  ///initialize edges
  Eigen::Matrix<typename DerivedV::Scalar, 3, 1> e[3];
  for (int k=0;k<3;k++)
    e[k] = V.row(v[(k+2)%3]) - V.row(v[(k+1)%3]);

  ///then consider area but also considering scale factor dur to overlaps

  double areaT = doublearea(f)/2.0;

  for (int x=0;x<3;x++)
    for (int y=0;y<3;y++)
      if (x!=y)
      {
        double num =  (e[x].dot(e[y]));
        val[x][y]  =  num/(4.0*areaT);
        val[x][y]  *= Handle_Stiffness[f];//f->stiffening;
      }

  ///set the matrix as diagonal
  SetDiagonal(val);
}

///initialize the RHS for a given face
///for minimization of Dirichlet's energy
template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::perElementRHS(int f,
                                                           double b[6],
                                                           double vector_field_scale)
{

  /// then set the rhs
  Eigen::Matrix<typename DerivedV::Scalar, 3, 1> scaled_Kreal;
  Eigen::Matrix<typename DerivedV::Scalar, 3, 1> scaled_Kimag;
  Eigen::Matrix<typename DerivedV::Scalar, 3, 1> fNorm = N.row(f);
  Eigen::Matrix<typename DerivedV::Scalar, 3, 1> p[3];
  p[0] = V.row(F(f,0));
  p[1] = V.row(F(f,1));
  p[2] = V.row(F(f,2));

  Eigen::Matrix<typename DerivedV::Scalar, 3, 1> neg_t[3];
  neg_t[0] = fNorm.cross(p[2] - p[1]);
  neg_t[1] = fNorm.cross(p[0] - p[2]);
  neg_t[2] = fNorm.cross(p[1] - p[0]);

  Eigen::Matrix<typename DerivedV::Scalar, 3, 1> K1,K2;
  K1 = PD1.row(f);
  K2 = -PD2.row(f); // TODO: the "-" accounts for the orientation of local_basis.h, adapt the code before and remove the "-"

  scaled_Kreal = K1*(vector_field_scale)/2;
  scaled_Kimag = K2*(vector_field_scale)/2;

  double stiff_val = Handle_Stiffness[f];

  b[0] = scaled_Kreal.dot(neg_t[0]) * stiff_val;
  b[1] = scaled_Kimag.dot(neg_t[0]) * stiff_val;
  b[2] = scaled_Kreal.dot(neg_t[1]) * stiff_val;
  b[3] = scaled_Kimag.dot(neg_t[1]) * stiff_val;
  b[4] = scaled_Kreal.dot(neg_t[2]) * stiff_val;
  b[5] = scaled_Kimag.dot(neg_t[2]) * stiff_val;

  //    if (f == 0)
  //    {
  //      cerr << "DEBUG!!!" << endl;
  //
  //
  //      for (unsigned z = 0; z<6; ++z)
  //        cerr << b[z] << " ";
  //      cerr << endl;
  //
  //      scaled_Kreal = K1*(vector_field_scale)/2;
  //      scaled_Kimag = -K2*(vector_field_scale)/2;
  //
  //      double stiff_val = Handle_Stiffness[f];
  //
  //      b[0] = scaled_Kreal.dot(neg_t[0]) * stiff_val;
  //      b[1] = scaled_Kimag.dot(neg_t[0]) * stiff_val;
  //      b[2] = scaled_Kreal.dot(neg_t[1]) * stiff_val;
  //      b[3] = scaled_Kimag.dot(neg_t[1]) * stiff_val;
  //      b[4] = scaled_Kreal.dot(neg_t[2]) * stiff_val;
  //      b[5] = scaled_Kimag.dot(neg_t[2]) * stiff_val;
  //
  //      for (unsigned z = 0; z<6; ++z)
  //        cerr << b[z] << " ";
  //      cerr << endl;
  //
  //    }

}

///evaluate the LHS and RHS for a single face
///for minimization of Dirichlet's energy
template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::PerElementSystemReal(int f,
                                                                  double val[3][3],
                                                                  int index[3][3][2],
                                                                  double b[6],
                                                                  double vector_field_scale)
{
  perElementLHS(f,val,index);
  perElementRHS(f,b,vector_field_scale);
}
///END ENERGY MINIMIZATION PART

///START FIXING VERTICES
///set a given vertex as fixed
template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::AddFixedVertex(int v)
{
  n_fixed_vars++;
  Hard_constraints.push_back(v);
}

///find vertex to fix in case we're using
///a vector field NB: multiple components not handled
template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::FindFixedVertField()
{
  Hard_constraints.clear();

  n_fixed_vars=0;
  ///fix the first singularity
  for (unsigned int v=0;v<V.rows();v++)
  {
    if (Handle_Singular(v))
    {
      AddFixedVertex(v);
      UV.row(v) << 0,0;
      return;
    }
  }

  ///if anything fixed fix the first
  AddFixedVertex(0); // TODO HERE IT ISSSSSS
  UV.row(0) << 0,0;
  std::cerr << "No vertices to fix, I am fixing the first vertex to the origin!" << std::endl;
}

///find hard constraint depending if using or not
///a vector field
template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::FindFixedVert()
{
  Hard_constraints.clear();
  FindFixedVertField();
}

template <typename DerivedV, typename DerivedF>
IGL_INLINE int igl::PoissonSolver<DerivedV, DerivedF>::GetFirstVertexIndex(int v)
{
  return HandleS_Index(VF[v][0],VFi[v][0]);
}

///fix the vertices which are flagged as fixed
template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::FixBlockedVertex()
{
  int offset_row = n_vert_vars*2 + num_cut_constraint*2;

  unsigned int constr_num = 0;
  for (unsigned int i=0;i<Hard_constraints.size();i++)
  {
    int v = Hard_constraints[i];

    ///get first index of the vertex that must blocked
    //int index=v->vertex_index[0];
    int index = GetFirstVertexIndex(v);

    ///multiply times 2 because of uv
    int indexvert = index*2;

    ///find the first free row to add the constraint
    int indexRow = (offset_row+constr_num*2);
    int indexCol = indexRow;

    ///add fixing constraint LHS
    AddValA(indexRow,indexvert,1);
    AddValA(indexRow+1,indexvert+1,1);

    ///add fixing constraint RHS
    AddValB(indexCol,  UV(v,0));
    AddValB(indexCol+1,UV(v,1));

    constr_num++;
  }
  assert(constr_num==n_fixed_vars);
}
///END FIXING VERTICES

///HANDLING SINGULARITY
//set the singularity round to integer location
template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::AddSingularityRound()
{
  for (unsigned int v=0;v<V.rows();v++)
  {
    if (Handle_Singular(v))
    {
      int index0=GetFirstVertexIndex(v);
      ids_to_round.push_back( index0*2   );
      ids_to_round.push_back((index0*2)+1);
    }
  }
}

template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::AddToRoundVertices(std::vector<int> ids)
{
  for (size_t i = 0; i < ids.size(); ++i)
  {
    if (ids[i] < 0 || ids[i] >= V.rows())
      std::cerr << "WARNING: Ignored round vertex constraint, vertex " << ids[i] << " does not exist in the mesh." << std::endl;
    int index0 = GetFirstVertexIndex(ids[i]);
    ids_to_round.push_back( index0*2   );
    ids_to_round.push_back((index0*2)+1);
  }
}

///START GENERIC SYSTEM FUNCTIONS
//build the laplacian matrix cyclyng over all rangemaps
//and over all faces
template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::BuildLaplacianMatrix(double vfscale)
{
  ///then for each face
  for (unsigned int f=0;f<F.rows();f++)
  {

    int var_idx[3]; //vertex variable indices

    for(int k = 0; k < 3; ++k)
      var_idx[k] = HandleS_Index(f,k);

    ///block of variables
    double val[3][3];
    ///block of vertex indexes
    int index[3][3][2];
    ///righe hand side
    double b[6];
    ///compute the system for the given face
    PerElementSystemReal(f, val,index, b, vfscale);

    //Add the element to the matrix
    Add33Block(val,index);

    ///add right hand side
    AddRHS(b,var_idx);
  }
}

///find different sized of the system
template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::FindSizes()
{
  ///find the vertex that need to be fixed
  FindFixedVert();

  ///REAL PART
  n_vert_vars = Handle_SystemInfo.num_vert_variables;

  ///INTEGER PART
  ///the total number of integer variables
  n_integer_vars = Handle_SystemInfo.num_integer_cuts;

  ///CONSTRAINT PART
  num_cut_constraint = Handle_SystemInfo.EdgeSeamInfo.size()*2;

  num_constraint_equations = num_cut_constraint*2 + n_fixed_vars*2 + num_userdefined_constraint;

  ///total variable of the system
  num_total_vars = n_vert_vars*2+n_integer_vars*2;

  ///initialize matrix size

  system_size = num_total_vars + num_constraint_equations;

  if (DEBUGPRINT)     printf("\n*** SYSTEM VARIABLES *** \n");
  if (DEBUGPRINT)     printf("* NUM REAL VERTEX VARIABLES %d \n",n_vert_vars);

  if (DEBUGPRINT)     printf("\n*** SINGULARITY *** \n ");
  if (DEBUGPRINT)     printf("* NUM SINGULARITY %d\n",(int)ids_to_round.size()/2);

  if (DEBUGPRINT)     printf("\n*** INTEGER VARIABLES *** \n");
  if (DEBUGPRINT)     printf("* NUM INTEGER VARIABLES %d \n",(int)n_integer_vars);

  if (DEBUGPRINT)     printf("\n*** CONSTRAINTS *** \n ");
  if (DEBUGPRINT)     printf("* NUM FIXED CONSTRAINTS %d\n",n_fixed_vars);
  if (DEBUGPRINT)     printf("* NUM CUTS CONSTRAINTS %d\n",num_cut_constraint);
  if (DEBUGPRINT)     printf("* NUM USER DEFINED CONSTRAINTS %d\n",num_userdefined_constraint);

  if (DEBUGPRINT)     printf("\n*** TOTAL SIZE *** \n");
  if (DEBUGPRINT)     printf("* TOTAL VARIABLE SIZE (WITH INTEGER TRASL) %d \n",num_total_vars);
  if (DEBUGPRINT)     printf("* TOTAL CONSTRAINTS %d \n",num_constraint_equations);
  if (DEBUGPRINT)     printf("* MATRIX SIZE  %d \n",system_size);
}

template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::AllocateSystem()
{
  S.initialize(system_size, system_size);
  printf("\n INITIALIZED SPARSE MATRIX OF %d x %d \n",system_size, system_size);
}

///intitialize the whole matrix
template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::InitMatrix()
{
  FindSizes();
  AllocateSystem();
}

///map back coordinates after that
///the system has been solved
template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::MapCoords()
{
  ///map coords to faces
  for (unsigned int f=0;f<F.rows();f++)
  {

    for (int k=0;k<3;k++)
    {
      //get the index of the variable in the system
      int indexUV = HandleS_Index(f,k);
      ///then get U and V coords
      double U=X[indexUV*2];
      double V=X[indexUV*2+1];

      WUV(f,k*2 + 0) = U;
      WUV(f,k*2 + 1) = V;
    }
  }

#if 0
  ///initialize the vector of integer variables to return their values
  Handle_SystemInfo.IntegerValues.resize(n_integer_vars*2);
  int baseIndex = (n_vert_vars)*2;
  int endIndex  = baseIndex+n_integer_vars*2;
  int index     = 0;
  for (int i=baseIndex; i<endIndex; i++)
  {
    ///assert that the value is an integer value
    double value=X[i];
    double diff = value-(int)floor(value+0.5);
    assert(diff<0.00000001);
    Handle_SystemInfo.IntegerValues[index] = value;
    index++;
  }
#endif
}

///END GENERIC SYSTEM FUNCTIONS

///set the constraints for the inter-range cuts
template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::BuildSeamConstraintsExplicitTranslation()
{
  ///add constraint(s) for every seam edge (not halfedge)
  int offset_row = n_vert_vars;
  ///current constraint row
  int constr_row = offset_row;
  ///current constraint
  unsigned int constr_num = 0;

  for (unsigned int i=0; i<num_cut_constraint/2; i++)
  {
    unsigned char interval = Handle_SystemInfo.EdgeSeamInfo[i].MMatch;
    if (interval==1)
      interval=3;
    else
      if(interval==3)
        interval=1;

    int p0  = Handle_SystemInfo.EdgeSeamInfo[i].v0;
    int p1  = Handle_SystemInfo.EdgeSeamInfo[i].v1;
    int p0p = Handle_SystemInfo.EdgeSeamInfo[i].v0p;
    int p1p = Handle_SystemInfo.EdgeSeamInfo[i].v1p;

    std::complex<double> rot = GetRotationComplex(interval);

    ///get the integer variable
    int integerVar = offset_row+Handle_SystemInfo.EdgeSeamInfo[i].integerVar;

    if (integer_rounding)
    {
      ids_to_round.push_back(integerVar*2);
      ids_to_round.push_back(integerVar*2+1);
    }

    AddComplexA(constr_row, p0 ,  rot);
    AddComplexA(constr_row, p0p,   -1);
    ///then translation...considering the rotation
    ///due to substitution
    AddComplexA(constr_row, integerVar, 1);

    AddValB(2*constr_row  ,0);
    AddValB(2*constr_row+1,0);
    constr_row +=1;
    constr_num++;

    AddComplexA(constr_row, p1,  rot);
    AddComplexA(constr_row, p1p, -1);

    ///other translation
    AddComplexA(constr_row, integerVar  , 1);

    AddValB(2*constr_row,0);
    AddValB(2*constr_row+1,0);

    constr_row +=1;
    constr_num++;
  }
}

///set the constraints for the inter-range cuts
template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::BuildUserDefinedConstraints()
{
  /// the user defined constraints are at the end
  int offset_row = n_vert_vars*2 + num_cut_constraint*2 + n_fixed_vars*2;

  ///current constraint row
  int constr_row = offset_row;

  assert(num_userdefined_constraint == userdefined_constraints.size());

  for (unsigned int i=0; i<num_userdefined_constraint; i++)
  {
    for (unsigned int j=0; j<userdefined_constraints[i].size()-1; ++j)
    {
      AddValA(constr_row, j ,  userdefined_constraints[i][j]);
    }

    AddValB(constr_row,userdefined_constraints[i][userdefined_constraints[i].size()-1]);

    constr_row +=1;
  }
}

///call of the mixed integer solver
template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::MixedIntegerSolve(double cone_grid_res,
                                                               bool direct_round,
                                                               int localIter)
{
  X = std::vector<double>((n_vert_vars+n_integer_vars)*2);

  ///variables part
  int ScalarSize = n_vert_vars*2;
  int SizeMatrix = (n_vert_vars+n_integer_vars)*2;

  if (DEBUGPRINT)
    printf("\n ALLOCATED X \n");

  ///matrix A
  gmm::col_matrix< gmm::wsvector< double > > A(SizeMatrix,SizeMatrix); // lhs matrix variables +

  ///constraints part
  int CsizeX = num_constraint_equations;
  int CsizeY = SizeMatrix+1;
  gmm::row_matrix< gmm::wsvector< double > > C(CsizeX,CsizeY); // constraints

  if (DEBUGPRINT)
    printf("\n ALLOCATED QMM STRUCTURES \n");

  std::vector<double> rhs(SizeMatrix,0);  // rhs

  if (DEBUGPRINT)
    printf("\n ALLOCATED RHS STRUCTURES \n");

  //// copy LHS
  for(int i = 0; i < (int)S.A().nentries(); ++i)
  {
    int row  = S.A().rowind()[i];
    int col  = S.A().colind()[i];
    int size =(int)S.nrows();
    assert(0 <= row && row < size);
    assert(0 <= col && col < size);

    // it's either part of the matrix
    if (row < ScalarSize)
    {
      A(row, col) += S.A().vals()[i];
    }
    // or it's a part of the constraint
    else
    {
      assert ((unsigned int)row < (n_vert_vars+num_constraint_equations)*2);
      int r = row - ScalarSize;
      assert(r   < CsizeX);
      assert(col < CsizeY);
      C(r  , col  ) +=  S.A().vals()[i];
    }
  }

  if (DEBUGPRINT)
    printf("\n SET %d INTEGER VALUES \n",n_integer_vars);

  ///add penalization term for integer variables
  double penalization = 0.000001;
  int offline_index   = ScalarSize;
  for(unsigned int i = 0; i < (n_integer_vars)*2; ++i)
  {
    int index=offline_index+i;
    A(index,index)=penalization;
  }

  if (DEBUGPRINT)
    printf("\n SET RHS \n");

  // copy RHS
  for(int i = 0; i < (int)ScalarSize; ++i)
  {
    rhs[i] = S.getRHSReal(i) * cone_grid_res;
  }

  // copy constraint RHS
  if (DEBUGPRINT)
    printf("\n SET %d CONSTRAINTS \n",num_constraint_equations);

  for(unsigned int i = 0; i < num_constraint_equations; ++i)
  {
    C(i, SizeMatrix) = -S.getRHSReal(ScalarSize + i) * cone_grid_res;
  }

  ///copy values back into S
  COMISO::ConstrainedSolver solver;

  solver.misolver().set_local_iters(localIter);
  solver.misolver().set_direct_rounding(direct_round);

  std::sort(ids_to_round.begin(),ids_to_round.end());
  std::vector<int>::iterator new_end=std::unique(ids_to_round.begin(),ids_to_round.end());
  int dist=distance(ids_to_round.begin(),new_end);
  ids_to_round.resize(dist);

  solver.solve( C, A, X, rhs, ids_to_round, 0.0, false, false);
}

template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::clearUserConstraint()
{
  num_userdefined_constraint = 0;
  userdefined_constraints.clear();
}

template <typename DerivedV, typename DerivedF>
IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::addSharpEdgeConstraint(int fid, int vid)
{
  // prepare constraint
  std::vector<int> c(Handle_SystemInfo.num_vert_variables*2 + 1);

  for (size_t i = 0; i < c.size(); ++i)
  {
    c[i] = 0;
  }

  int v1 = F(fid,vid);
  int v2 = F(fid,(vid+1)%3);

  Eigen::Matrix<typename DerivedV::Scalar, 3, 1> e = V.row(v2) - V.row(v1);

  int v1i = GetFirstVertexIndex(v1);
  int v2i = GetFirstVertexIndex(v2);

  double d1 = fabs(e.dot(PD1.row(fid)));
  double d2 = fabs(e.dot(PD2.row(fid)));

  int offset = 0;

  if (d1>d2)
    offset = 1;

  ids_to_round.push_back((v1i * 2) + offset);
  ids_to_round.push_back((v2i * 2) + offset);

  // add constraint
  c[(v1i * 2) + offset] =  1;
  c[(v2i * 2) + offset] = -1;

  // add to the user-defined constraints
  num_userdefined_constraint++;
  userdefined_constraints.push_back(c);

}



template <typename DerivedV, typename DerivedF, typename DerivedU>
IGL_INLINE igl::MIQ_class<DerivedV, DerivedF, DerivedU>::MIQ_class(const Eigen::PlainObjectBase<DerivedV> &V_,
                                                        const Eigen::PlainObjectBase<DerivedF> &F_,
                                                        const Eigen::PlainObjectBase<DerivedV> &PD1_combed,
                                                        const Eigen::PlainObjectBase<DerivedV> &PD2_combed,
                                                        const Eigen::PlainObjectBase<DerivedV> &BIS1_combed,
                                                        const Eigen::PlainObjectBase<DerivedV> &BIS2_combed,
                                                        const Eigen::Matrix<int, Eigen::Dynamic, 3> &Handle_MMatch,
                                                        const Eigen::Matrix<int, Eigen::Dynamic, 1> &Handle_Singular,
                                                        const Eigen::Matrix<int, Eigen::Dynamic, 1> &Handle_SingularDegree,
                                                        const Eigen::Matrix<int, Eigen::Dynamic, 3> &Handle_Seams,
                                                        Eigen::PlainObjectBase<DerivedU> &UV,
                                                        Eigen::PlainObjectBase<DerivedF> &FUV,
                                                        double GradientSize,
                                                        double Stiffness,
                                                        bool DirectRound,
                                                        int iter,
                                                        int localIter,
                                                        bool DoRound,
                                                        std::vector<int> roundVertices,
                                                        std::vector<std::vector<int> > hardFeatures):
V(V_),
F(F_)
{
  igl::local_basis(V,F,B1,B2,B3);
  igl::triangle_triangle_adjacency(V,F,TT,TTi);

  // Prepare indexing for the linear system
  VertexIndexing<DerivedV, DerivedF> VInd(V, F, TT, TTi, BIS1_combed, BIS2_combed, Handle_MMatch, Handle_Singular, Handle_SingularDegree, Handle_Seams);

  VInd.InitMapping();
  VInd.InitFaceIntegerVal();
  VInd.InitSeamInfo();

  Eigen::PlainObjectBase<DerivedV> PD1_combed_for_poisson, PD2_combed_for_poisson;
  // Rotate by 90 degrees CCW
  PD1_combed_for_poisson.setZero(BIS1_combed.rows(),3);
  PD2_combed_for_poisson.setZero(BIS2_combed.rows(),3);
  for (unsigned i=0; i<PD1_combed.rows();++i)
  {
    double n1 = PD1_combed.row(i).norm();
    double n2 = PD2_combed.row(i).norm();

    double a1 = atan2(B2.row(i).dot(PD1_combed.row(i)),B1.row(i).dot(PD1_combed.row(i)));
    double a2 = atan2(B2.row(i).dot(PD2_combed.row(i)),B1.row(i).dot(PD2_combed.row(i)));

    a1 += M_PI/2;
    a2 += M_PI/2;


    PD1_combed_for_poisson.row(i) = cos(a1) * B1.row(i) + sin(a1) * B2.row(i);
    PD2_combed_for_poisson.row(i) = cos(a2) * B1.row(i) + sin(a2) * B2.row(i);

    PD1_combed_for_poisson.row(i) = PD1_combed_for_poisson.row(i).normalized() * n1;
    PD2_combed_for_poisson.row(i) = PD2_combed_for_poisson.row(i).normalized() * n2;
  }


  // Assemble the system and solve
  PoissonSolver<DerivedV, DerivedF> PSolver(V,
                                            F,
                                            TT,
                                            TTi,
                                            PD1_combed_for_poisson,
                                            PD2_combed_for_poisson,
                                            VInd.HandleS_Index,
                                            VInd.Handle_Singular,
                                            VInd.Handle_SystemInfo);
  Handle_Stiffness = Eigen::VectorXd::Constant(F.rows(),1);


  if (iter > 0) // do stiffening
  {
    for (int i=0;i<iter;i++)
    {
      PSolver.SolvePoisson(Handle_Stiffness, GradientSize,1.f,DirectRound,localIter,DoRound,roundVertices,hardFeatures);
      int nflips=NumFlips(PSolver.WUV);
      bool folded = updateStiffeningJacobianDistorsion(GradientSize,PSolver.WUV);
      printf("ITERATION %d FLIPS %d \n",i,nflips);
      if (!folded)break;
    }
  }
  else
  {
    PSolver.SolvePoisson(Handle_Stiffness,GradientSize,1.f,DirectRound,localIter,DoRound,roundVertices,hardFeatures);
  }

  int nflips=NumFlips(PSolver.WUV);
  printf("**** END OPTIMIZING #FLIPS %d  ****\n",nflips);

  fflush(stdout);
  WUV = PSolver.WUV;

}

template <typename DerivedV, typename DerivedF, typename DerivedU>
IGL_INLINE void igl::MIQ_class<DerivedV, DerivedF, DerivedU>::extractUV(Eigen::PlainObjectBase<DerivedU> &UV_out,
                                                             Eigen::PlainObjectBase<DerivedF> &FUV_out)
{
  //      int f = F.rows();
  int f = WUV.rows();

  unsigned vtfaceid[f*3];
  std::vector<double> vtu;
  std::vector<double> vtv;

  std::vector<std::vector<double> > listUV;
  unsigned counter = 0;

  for (unsigned i=0; i<f; ++i)
  {
    for (unsigned j=0; j<3; ++j)
    {
      std::vector<double> t(3);
      t[0] = WUV(i,j*2 + 0);
      t[1] = WUV(i,j*2 + 1);
      t[2] = counter++;
      listUV.push_back(t);
    }
  }
  std::sort(listUV.begin(),listUV.end());

  counter = 0;
  unsigned k = 0;
  while (k < f*3)
  {
    double u = listUV[k][0];
    double v = listUV[k][1];
    unsigned id = round(listUV[k][2]);

    vtfaceid[id] = counter;
    vtu.push_back(u);
    vtv.push_back(v);

    unsigned j=1;
    while(k+j < f*3 && u == listUV[k+j][0] && v == listUV[k+j][1])
    {
      unsigned tid = round(listUV[k+j][2]);
      vtfaceid[tid] = counter;
      ++j;
    }
    k = k+j;
    counter++;
  }

  UV_out.resize(vtu.size(),2);
  for (unsigned i=0; i<vtu.size(); ++i)
  {
    UV_out(i,0) = vtu[i];
    UV_out(i,1) = vtv[i];
  }

  FUV_out.resize(f,3);

  unsigned vcounter = 0;
  for (unsigned i=0; i<f; ++i)
  {
    FUV_out(i,0)  = vtfaceid[vcounter++];
    FUV_out(i,1)  = vtfaceid[vcounter++];
    FUV_out(i,2)  = vtfaceid[vcounter++];
  }

}

template <typename DerivedV, typename DerivedF, typename DerivedU>
IGL_INLINE int igl::MIQ_class<DerivedV, DerivedF, DerivedU>::NumFlips(const Eigen::MatrixXd& WUV)
{
  int numFl=0;
  for (unsigned int i=0;i<F.rows();i++)
  {
    if (IsFlipped(i, WUV))
      numFl++;
  }
  return numFl;
}

template <typename DerivedV, typename DerivedF, typename DerivedU>
IGL_INLINE double igl::MIQ_class<DerivedV, DerivedF, DerivedU>::Distortion(int f, double h, const Eigen::MatrixXd& WUV)
{
  assert(h > 0);

  Eigen::Vector2d uv0,uv1,uv2;

  uv0 << WUV(f,0), WUV(f,1);
  uv1 << WUV(f,2), WUV(f,3);
  uv2 << WUV(f,4), WUV(f,5);

  Eigen::Matrix<typename DerivedV::Scalar, 3, 1> p0 = V.row(F(f,0));
  Eigen::Matrix<typename DerivedV::Scalar, 3, 1> p1 = V.row(F(f,1));
  Eigen::Matrix<typename DerivedV::Scalar, 3, 1> p2 = V.row(F(f,2));

  Eigen::Matrix<typename DerivedV::Scalar, 3, 1> norm = (p1 - p0).cross(p2 - p0);
  double area2 = norm.norm();
  double area2_inv  = 1.0 / area2;
  norm *= area2_inv;

  if (area2 > 0)
  {
    // Singular values of the Jacobian
    Eigen::Matrix<typename DerivedV::Scalar, 3, 1> neg_t0 = norm.cross(p2 - p1);
    Eigen::Matrix<typename DerivedV::Scalar, 3, 1> neg_t1 = norm.cross(p0 - p2);
    Eigen::Matrix<typename DerivedV::Scalar, 3, 1> neg_t2 = norm.cross(p1 - p0);

    Eigen::Matrix<typename DerivedV::Scalar, 3, 1> diffu =  (neg_t0 * uv0(0) +neg_t1 *uv1(0) +  neg_t2 * uv2(0) )*area2_inv;
    Eigen::Matrix<typename DerivedV::Scalar, 3, 1> diffv = (neg_t0 * uv0(1) + neg_t1*uv1(1) +  neg_t2*uv2(1) )*area2_inv;

    // first fundamental form
    double I00 = diffu.dot(diffu);  // guaranteed non-neg
    double I01 = diffu.dot(diffv);  // I01 = I10
    double I11 = diffv.dot(diffv);  // guaranteed non-neg

    // eigenvalues of a 2x2 matrix
    // [a00 a01]
    // [a10 a11]
    // 1/2 * [ (a00 + a11) +/- sqrt((a00 - a11)^2 + 4 a01 a10) ]
    double trI = I00 + I11;                     // guaranteed non-neg
    double diffDiag = I00 - I11;                // guaranteed non-neg
    double sqrtDet = sqrt(std::max(0.0, diffDiag*diffDiag +
                                   4 * I01 * I01)); // guaranteed non-neg
    double sig1 = 0.5 * (trI + sqrtDet); // higher singular value
    double sig2 = 0.5 * (trI - sqrtDet); // lower singular value

    // Avoid sig2 < 0 due to numerical error
    if (fabs(sig2) < 1.0e-8)
      sig2 = 0;

    assert(sig1 >= 0);
    assert(sig2 >= 0);

    if (sig2 < 0) {
      printf("Distortion will be NaN! sig1^2 is negative (%lg)\n",
             sig2);
    }

    // The singular values of the Jacobian are the sqrts of the
    // eigenvalues of the first fundamental form.
    sig1 = sqrt(sig1);
    sig2 = sqrt(sig2);

    // distortion
    double tao = IsFlipped(f,WUV) ? -1 : 1;
    double factor = tao / h;
    double lam = fabs(factor * sig1 - 1) + fabs(factor * sig2 - 1);
    return lam;
  }
  else {
    return 10; // something "large"
  }
}

////////////////////////////////////////////////////////////////////////////
// Approximate the distortion laplacian using a uniform laplacian on
//  the dual mesh:
//      ___________
//      \-1 / \-1 /
//       \ / 3 \ /
//        \-----/
//         \-1 /
//          \ /
//
//  @param[in]  f   facet on which to compute distortion laplacian
//  @param[in]  h   scaling factor applied to cross field
//  @return     distortion laplacian for f
///////////////////////////////////////////////////////////////////////////
template <typename DerivedV, typename DerivedF, typename DerivedU>
IGL_INLINE double igl::MIQ_class<DerivedV, DerivedF, DerivedU>::LaplaceDistortion(const int f, double h, const Eigen::MatrixXd& WUV)
{
  double mydist = Distortion(f, h, WUV);
  double lapl=0;
  for (int i=0;i<3;i++)
  {
    if (TT(f,i) != -1)
      lapl += (mydist - Distortion(TT(f,i), h, WUV));
  }
  return lapl;
}

template <typename DerivedV, typename DerivedF, typename DerivedU>
IGL_INLINE bool igl::MIQ_class<DerivedV, DerivedF, DerivedU>::updateStiffeningJacobianDistorsion(double grad_size, const Eigen::MatrixXd& WUV)
{
  bool flipped = NumFlips(WUV)>0;

  if (!flipped)
    return false;

  double maxL=0;
  double maxD=0;

  if (flipped)
  {
    const double c = 1.0;
    const double d = 5.0;

    for (unsigned int i = 0; i < F.rows(); ++i)
    {
      double dist=Distortion(i,grad_size,WUV);
      if (dist > maxD)
        maxD=dist;

      double absLap=fabs(LaplaceDistortion(i, grad_size,WUV));
      if (absLap > maxL)
        maxL = absLap;

      double stiffDelta = std::min(c * absLap, d);

      Handle_Stiffness[i]+=stiffDelta;
    }
  }
  printf("Maximum Distorsion %4.4f \n",maxD);
  printf("Maximum Laplacian %4.4f \n",maxL);
  return flipped;
}

template <typename DerivedV, typename DerivedF, typename DerivedU>
IGL_INLINE bool igl::MIQ_class<DerivedV, DerivedF, DerivedU>::IsFlipped(const Eigen::Vector2d &uv0,
                                                                    const Eigen::Vector2d &uv1,
                                                                    const Eigen::Vector2d &uv2)
{
  Eigen::Vector2d e0 = (uv1-uv0);
  Eigen::Vector2d e1 = (uv2-uv0);

  double Area = e0(0)*e1(1) - e0(1)*e1(0);
  return (Area<=0);
}

template <typename DerivedV, typename DerivedF, typename DerivedU>
IGL_INLINE bool igl::MIQ_class<DerivedV, DerivedF, DerivedU>::IsFlipped(const int i, const Eigen::MatrixXd& WUV)
{
  Eigen::Vector2d uv0,uv1,uv2;
  uv0 << WUV(i,0), WUV(i,1);
  uv1 << WUV(i,2), WUV(i,3);
  uv2 << WUV(i,4), WUV(i,5);

  return (IsFlipped(uv0,uv1,uv2));
}




template <typename DerivedV, typename DerivedF, typename DerivedU>
IGL_INLINE void igl::miq(const Eigen::PlainObjectBase<DerivedV> &V,
                         const Eigen::PlainObjectBase<DerivedF> &F,
                         const Eigen::PlainObjectBase<DerivedV> &PD1_combed,
                         const Eigen::PlainObjectBase<DerivedV> &PD2_combed,
                         const Eigen::PlainObjectBase<DerivedV> &BIS1_combed,
                         const Eigen::PlainObjectBase<DerivedV> &BIS2_combed,
                         const Eigen::Matrix<int, Eigen::Dynamic, 3> &Handle_MMatch,
                         const Eigen::Matrix<int, Eigen::Dynamic, 1> &Handle_Singular,
                         const Eigen::Matrix<int, Eigen::Dynamic, 1> &Handle_SingularDegree,
                         const Eigen::Matrix<int, Eigen::Dynamic, 3> &Handle_Seams,
                         Eigen::PlainObjectBase<DerivedU> &UV,
                         Eigen::PlainObjectBase<DerivedF> &FUV,
                         double GradientSize,
                         double Stiffness,
                         bool DirectRound,
                         int iter,
                         int localIter,
                         bool DoRound,
                         std::vector<int> roundVertices,
                         std::vector<std::vector<int> > hardFeatures)
{
  GradientSize = GradientSize/(V.colwise().maxCoeff()-V.colwise().minCoeff()).norm();

  igl::MIQ_class<DerivedV, DerivedF, DerivedU> miq(V,
                                                   F,
                                                   PD1_combed,
                                                   PD2_combed,
                                                   BIS1_combed,
                                                   BIS2_combed,
                                                   Handle_MMatch,
                                                   Handle_Singular,
                                                   Handle_SingularDegree,
                                                   Handle_Seams,
                                                   UV,
                                                   FUV,
                                                   GradientSize,
                                                   Stiffness,
                                                   DirectRound,
                                                   iter,
                                                   localIter,
                                                   DoRound,
                                                   roundVertices,
                                                   hardFeatures);

  miq.extractUV(UV,FUV);
}

template <typename DerivedV, typename DerivedF, typename DerivedU>
IGL_INLINE void igl::miq(const Eigen::PlainObjectBase<DerivedV> &V,
                         const Eigen::PlainObjectBase<DerivedF> &F,
                         const Eigen::PlainObjectBase<DerivedV> &PD1,
                         const Eigen::PlainObjectBase<DerivedV> &PD2,
                         Eigen::PlainObjectBase<DerivedU> &UV,
                         Eigen::PlainObjectBase<DerivedF> &FUV,
                         double GradientSize,
                         double Stiffness,
                         bool DirectRound,
                         int iter,
                         int localIter,
                         bool DoRound,
                         std::vector<int> roundVertices,
                         std::vector<std::vector<int> > hardFeatures)
{
  Eigen::MatrixXd PD2i = PD2;
  if (PD2i.size() == 0)
  {
    Eigen::MatrixXd B1, B2, B3;
    igl::local_basis(V,F,B1,B2,B3);
    PD2i = igl::rotate_vectors(V,Eigen::VectorXd::Constant(1,M_PI/2),B1,B2);
  }

  Eigen::PlainObjectBase<DerivedV> BIS1, BIS2;
  igl::compute_frame_field_bisectors(V, F, PD1, PD2, BIS1, BIS2);

  Eigen::PlainObjectBase<DerivedV> BIS1_combed, BIS2_combed;
  igl::comb_cross_field(V, F, BIS1, BIS2, BIS1_combed, BIS2_combed);

  Eigen::Matrix<int, Eigen::Dynamic, 3> Handle_MMatch;
  igl::cross_field_missmatch(V, F, BIS1_combed, BIS2_combed, true, Handle_MMatch);

  Eigen::Matrix<int, Eigen::Dynamic, 1> isSingularity, singularityIndex;
  igl::find_cross_field_singularities(V, F, Handle_MMatch, isSingularity, singularityIndex);

  Eigen::Matrix<int, Eigen::Dynamic, 3> Handle_Seams;
  igl::cut_mesh_from_singularities(V, F, Handle_MMatch, isSingularity, singularityIndex, Handle_Seams);

  Eigen::PlainObjectBase<DerivedV> PD1_combed, PD2_combed;
  igl::comb_frame_field(V, F, PD1, PD2, BIS1_combed, BIS2_combed, PD1_combed, PD2_combed);

  igl::miq(V,
           F,
           PD1_combed,
           PD2_combed,
           BIS1_combed,
           BIS2_combed,
           Handle_MMatch,
           isSingularity,
           singularityIndex,
           Handle_Seams,
           UV,
           FUV,
           GradientSize,
           Stiffness,
           DirectRound,
           iter,
           localIter,
           DoRound,
           roundVertices,
           hardFeatures);

}

#ifdef IGL_STATIC_LIBRARY
// Explicit template specialization
template void igl::mixed_integer_quadrangulate<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 2, 0, -1, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, double, double, bool, int, int, bool, std::__1::vector<int, std::__1::allocator<int> >, std::__1::vector<std::__1::vector<int, std::__1::allocator<int> >, std::__1::allocator<std::__1::vector<int, std::__1::allocator<int> > > >);
#endif