Fix CI Build: avoid test on windows, remove Comiso module, specific xcode on github (#2384)

* just avoid failing test on windows 😔

* mayfil instead

* mac os x xcode bug fix

* rm bonus endif

* rm comiso

.github/workflows/continuous.yml

@@ -52,6 +52,12 @@ jobs:
         run: |
           HOMEBREW_NO_AUTO_UPDATE=1 brew install ccache
 
+      - name: Setup Xcode version
+        if: runner.os == 'macOS'
+        uses: maxim-lobanov/setup-xcode@v1
+        with:
+          xcode-version: '15.2'
+
       - name: Cache Build
         id: cache-build
         uses: actions/cache@v2

cmake/igl/modules/copyleft/comiso.cmake

@@ -1,27 +0,0 @@
-# 1. Define module
-igl_add_library(igl_copyleft_comiso)
-
-# 2. Include headers
-include(GNUInstallDirs)
-target_include_directories(igl_copyleft_comiso ${IGL_SCOPE}
-    $<BUILD_INTERFACE:${libigl_SOURCE_DIR}/include>
-    $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>
-)
-
-# 3. Target sources
-file(GLOB INC_FILES "${libigl_SOURCE_DIR}/include/igl/copyleft/comiso/*.h")
-file(GLOB SRC_FILES "${libigl_SOURCE_DIR}/include/igl/copyleft/comiso/*.cpp")
-igl_target_sources(igl_copyleft_comiso ${INC_FILES} ${SRC_FILES})
-
-# 4. Dependencies
-include(comiso)
-igl_include(copyleft core)
-target_link_libraries(igl_copyleft_comiso ${IGL_SCOPE}
-    igl::core
-    igl_copyleft::core
-    CoMISo::CoMISo
-)
-
-# 5. Unit tests
-file(GLOB SRC_FILES "${libigl_SOURCE_DIR}/tests/include/igl/copyleft/comiso/*.cpp")
-igl_add_test(igl_copyleft_comiso ${SRC_FILES})

cmake/libigl.cmake

@@ -25,7 +25,6 @@ igl_include_optional(xml)
 # Libigl copyleft modules
 igl_include_optional(copyleft core)
 igl_include_optional(copyleft cgal)
-igl_include_optional(copyleft comiso)
 igl_include_optional(copyleft tetgen)
 
 # Libigl restricted modules

cmake/recipes/external/comiso.cmake

@@ -1,31 +0,0 @@
-if(TARGET CoMISo::CoMISo)
-    return()
-endif()
-
-message(STATUS "Third-party: creating target 'CoMISo::CoMISo'")
-
-include(FetchContent)
-FetchContent_Declare(
-    comiso
-    GIT_REPOSITORY https://github.com/libigl/CoMISo.git
-    GIT_TAG 536440e714f412e7ef6c0b96b90ba37b1531bb39
-)
-
-include(eigen)
-
-FetchContent_MakeAvailable(comiso)
-
-add_library(CoMISo::CoMISo ALIAS CoMISo)
-
-# Copy .hh headers into a subfolder `CoMISo/`
-file(GLOB_RECURSE INC_FILES "${comiso_SOURCE_DIR}/*.hh" "${comiso_SOURCE_DIR}/*.cc")
-set(output_folder "${CMAKE_CURRENT_BINARY_DIR}/CoMISo/include/CoMISo")
-message(VERBOSE "Copying CoMISo headers to '${output_folder}'")
-foreach(filepath IN ITEMS ${INC_FILES})
-    file(RELATIVE_PATH filename "${comiso_SOURCE_DIR}" ${filepath})
-    configure_file(${filepath} "${output_folder}/${filename}" COPYONLY)
-endforeach()
-
-target_include_directories(CoMISo PUBLIC ${CMAKE_CURRENT_BINARY_DIR}/CoMISo/include)
-
-set_target_properties(CoMISo PROPERTIES FOLDER ThirdParty)

include/igl/copyleft/comiso/frame_field.cpp

@@ -1,688 +0,0 @@
-// This file is part of libigl, a simple c++ geometry processing library.
-//
-// Copyright (C) 2015 Daniele Panozzo <[email protected]>
-//
-// This Source Code Form is subject to the terms of the Mozilla Public License
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can
-// obtain one at http://mozilla.org/MPL/2.0/.
-#include "frame_field.h"
-
-#include "../../triangle_triangle_adjacency.h"
-#include "../../edge_topology.h"
-#include "../../per_face_normals.h"
-#include "nrosy.h"
-#include <iostream>
-
-namespace igl
-{
-namespace copyleft
-{
-namespace comiso
-{
-
-class FrameInterpolator
-{
-public:
-  // Init
-  IGL_INLINE FrameInterpolator(const Eigen::MatrixXd& _V, const Eigen::MatrixXi& _F);
-  IGL_INLINE ~FrameInterpolator();
-
-  // Reset constraints (at least one constraint must be present or solve will fail)
-  IGL_INLINE void resetConstraints();
-
-  IGL_INLINE void setConstraint(const int fid, const Eigen::VectorXd& v);
-
-  IGL_INLINE void interpolateSymmetric();
-
-  // Generate the frame field
-  IGL_INLINE void solve();
-
-  // Convert the frame field in the canonical representation
-  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(const Eigen::MatrixXd& TP, const double theta, const Eigen::VectorXd& S, Eigen::RowVectorXd& v);
-
-  IGL_INLINE Eigen::MatrixXd getFieldPerFace();
-
-  IGL_INLINE void PolarDecomposition(Eigen::MatrixXd V, Eigen::MatrixXd& U, Eigen::MatrixXd& P);
-
-  // Symmetric
-  Eigen::MatrixXd S;
-  std::vector<bool> S_c;
-
-  // -------------------------------------------------
-
-  // Face Topology
-  Eigen::MatrixXi TT, TTi;
-
-  // Two faces are consistent if their representative vector are taken modulo PI
-  std::vector<bool> edge_consistency;
-  Eigen::MatrixXi   edge_consistency_TT;
-
-private:
-  IGL_INLINE double mod2pi(double d);
-  IGL_INLINE double modpi2(double d);
-  IGL_INLINE double modpi(double d);
-
-  // Convert a direction on the tangent space into an angle
-  IGL_INLINE double vector2theta(const Eigen::MatrixXd& TP, const Eigen::RowVectorXd& v);
-
-  // Convert an angle in a vector in the tangent space
-  IGL_INLINE Eigen::RowVectorXd theta2vector(const Eigen::MatrixXd& TP, const double theta);
-
-  // Interpolate the cross field (theta)
-  IGL_INLINE void interpolateCross();
-
-  // Compute difference between reference frames
-  IGL_INLINE void computek();
-
-  // Compute edge consistency
-  IGL_INLINE void compute_edge_consistency();
-
-  // Cross field direction
-  Eigen::VectorXd thetas;
-  std::vector<bool> thetas_c;
-
-  // Edge Topology
-  Eigen::MatrixXi EV, FE, EF;
-  std::vector<bool> isBorderEdge;
-
-  // Angle between two reference frames
-  // R(k) * t0 = t1
-  Eigen::VectorXd k;
-
-  // Mesh
-  Eigen::MatrixXd V;
-  Eigen::MatrixXi F;
-
-  // Normals per face
-  Eigen::MatrixXd N;
-
-  // Reference frame per triangle
-  std::vector<Eigen::MatrixXd> TPs;
-
-};
-
-FrameInterpolator::FrameInterpolator(const Eigen::MatrixXd& _V, const Eigen::MatrixXi& _F)
-{
-  using namespace std;
-  using namespace Eigen;
-
-  V = _V;
-  F = _F;
-
-  assert(V.rows() > 0);
-  assert(F.rows() > 0);
-
-
-  // Generate topological relations
-  igl::triangle_triangle_adjacency(F,TT,TTi);
-  igl::edge_topology(V,F, EV, FE, EF);
-
-  // Flag border edges
-  isBorderEdge.resize(EV.rows());
-  for(unsigned i=0; i<EV.rows(); ++i)
-    isBorderEdge[i] = (EF(i,0) == -1) || ((EF(i,1) == -1));
-
-  // Generate normals per face
-  igl::per_face_normals(V, F, N);
-
-  // Generate reference frames
-  for(unsigned fid=0; fid<F.rows(); ++fid)
-  {
-    // First edge
-    Vector3d e1 = V.row(F(fid,1)) - V.row(F(fid,0));
-    e1.normalize();
-    Vector3d e2 = N.row(fid);
-    e2 = e2.cross(e1);
-    e2.normalize();
-
-    MatrixXd TP(2,3);
-    TP << e1.transpose(), e2.transpose();
-    TPs.push_back(TP);
-  }
-
-  // Reset the constraints
-  resetConstraints();
-
-  // Compute k, differences between reference frames
-  computek();
-
-  // Alloc internal variables
-  thetas            = VectorXd::Zero(F.rows());
-  S = MatrixXd::Zero(F.rows(),3);
-
-  compute_edge_consistency();
-}
-
-FrameInterpolator::~FrameInterpolator()
-{
-
-}
-
-double FrameInterpolator::mod2pi(double d)
-{
-  while(d<0)
-    d = d + (2.0*igl::PI);
-
-  return fmod(d, (2.0*igl::PI));
-}
-
-double FrameInterpolator::modpi2(double d)
-{
-  while(d<0)
-    d = d + (igl::PI/2.0);
-
-  return fmod(d, (igl::PI/2.0));
-}
-
-double FrameInterpolator::modpi(double d)
-{
-  while(d<0)
-    d = d + (igl::PI);
-
-  return fmod(d, (igl::PI));
-}
-
-
-double FrameInterpolator::vector2theta(const Eigen::MatrixXd& TP, const Eigen::RowVectorXd& v)
-{
-  // Project onto the tangent plane
-  Eigen::Vector2d vp = TP * v.transpose();
-
-  // Convert to angle
-  double theta = atan2(vp(1),vp(0));
-  return theta;
-}
-
-Eigen::RowVectorXd FrameInterpolator::theta2vector(const Eigen::MatrixXd& TP, const double theta)
-{
-  Eigen::Vector2d vp(cos(theta),sin(theta));
-  return vp.transpose() * TP;
-}
-
-void FrameInterpolator::interpolateCross()
-{
-  using namespace std;
-  using namespace Eigen;
-
-  //olga: was
-  // NRosyField nrosy(V,F);
-  // for (unsigned i=0; i<F.rows(); ++i)
-    // if(thetas_c[i])
-      // nrosy.setConstraintHard(i,theta2vector(TPs[i],thetas(i)));
-  // nrosy.solve(4);
-  // MatrixXd R = nrosy.getFieldPerFace();
-
-  //olga: is
-  Eigen::MatrixXd R;
-  Eigen::VectorXd S;
-  Eigen::VectorXi b; b.resize(F.rows(),1);
-  Eigen::MatrixXd bc; bc.resize(F.rows(),3);
-  int num = 0;
-  for (unsigned i=0; i<F.rows(); ++i)
-    if(thetas_c[i])
-      {
-        b[num] = i;
-        bc.row(num) = theta2vector(TPs[i],thetas(i));
-        num++;
-      }
-  b.conservativeResize(num,Eigen::NoChange);
-  bc.conservativeResize(num,Eigen::NoChange);
-
-  igl::copyleft::comiso::nrosy(V, F, b, bc, 4, R, S);
-  //olga:end
-  assert(R.rows() == F.rows());
-
-  for (unsigned i=0; i<F.rows(); ++i)
-    thetas(i) = vector2theta(TPs[i],R.row(i));
-}
-
-void FrameInterpolator::resetConstraints()
-{
-  thetas_c.resize(F.rows());
-  S_c.resize(F.rows());
-
-  for(unsigned i=0; i<F.rows(); ++i)
-  {
-    thetas_c[i]  = false;
-    S_c[i] = false;
-  }
-
-}
-
-void FrameInterpolator::compute_edge_consistency()
-{
-  using namespace std;
-  using namespace Eigen;
-
-  // Compute per-edge consistency
-  edge_consistency.resize(EF.rows());
-  edge_consistency_TT = MatrixXi::Constant(TT.rows(),3,-1);
-
-  // For every non-border edge
-  for (unsigned eid=0; eid<EF.rows(); ++eid)
-  {
-    if (!isBorderEdge[eid])
-    {
-      int fid0 = EF(eid,0);
-      int fid1 = EF(eid,1);
-
-      double theta0 = thetas(fid0);
-      double theta1 = thetas(fid1);
-
-      theta0 = theta0 + k(eid);
-
-      double r = modpi(theta0-theta1);
-
-      edge_consistency[eid] = r < igl::PI/4.0 || r > 3*(igl::PI/4.0);
-
-      // Copy it into edge_consistency_TT
-      int i1 = -1;
-      int i2 = -1;
-      for (unsigned i=0; i<3; ++i)
-      {
-        if (TT(fid0,i) == fid1)
-          i1 = i;
-        if (TT(fid1,i) == fid0)
-          i2 = i;
-      }
-      assert(i1 != -1);
-      assert(i2 != -1);
-
-      edge_consistency_TT(fid0,i1) = edge_consistency[eid];
-      edge_consistency_TT(fid1,i2) = edge_consistency[eid];
-    }
-  }
-}
-
-void FrameInterpolator::computek()
-{
-  using namespace std;
-  using namespace Eigen;
-
-  k.resize(EF.rows());
-
-  // For every non-border edge
-  for (unsigned eid=0; eid<EF.rows(); ++eid)
-  {
-    if (!isBorderEdge[eid])
-    {
-      int fid0 = EF(eid,0);
-      int fid1 = EF(eid,1);
-
-      Vector3d N0 = N.row(fid0);
-      //Vector3d N1 = N.row(fid1);
-
-      // find common edge on triangle 0 and 1
-      int fid0_vc = -1;
-      int fid1_vc = -1;
-      for (unsigned i=0;i<3;++i)
-      {
-        if (EV(eid,0) == F(fid0,i))
-          fid0_vc = i;
-        if (EV(eid,1) == F(fid1,i))
-          fid1_vc = i;
-      }
-      assert(fid0_vc != -1);
-      assert(fid1_vc != -1);
-
-      Vector3d common_edge = V.row(F(fid0,(fid0_vc+1)%3)) - V.row(F(fid0,fid0_vc));
-      common_edge.normalize();
-
-      // Map the two triangles in a new space where the common edge is the x axis and the N0 the z axis
-      MatrixXd P(3,3);
-      VectorXd o = V.row(F(fid0,fid0_vc));
-      VectorXd tmp = -N0.cross(common_edge);
-      P << common_edge, tmp, N0;
-      P.transposeInPlace();
-
-
-      MatrixXd V0(3,3);
-      V0.row(0) = V.row(F(fid0,0)).transpose() -o;
-      V0.row(1) = V.row(F(fid0,1)).transpose() -o;
-      V0.row(2) = V.row(F(fid0,2)).transpose() -o;
-
-      V0 = (P*V0.transpose()).transpose();
-
-      assert(V0(0,2) < 10e-10);
-      assert(V0(1,2) < 10e-10);
-      assert(V0(2,2) < 10e-10);
-
-      MatrixXd V1(3,3);
-      V1.row(0) = V.row(F(fid1,0)).transpose() -o;
-      V1.row(1) = V.row(F(fid1,1)).transpose() -o;
-      V1.row(2) = V.row(F(fid1,2)).transpose() -o;
-      V1 = (P*V1.transpose()).transpose();
-
-      assert(V1(fid1_vc,2) < 10e-10);
-      assert(V1((fid1_vc+1)%3,2) < 10e-10);
-
-      // compute rotation R such that R * N1 = N0
-      // i.e. map both triangles to the same plane
-      double alpha = -atan2(V1((fid1_vc+2)%3,2),V1((fid1_vc+2)%3,1));
-
-      MatrixXd R(3,3);
-      R << 1,          0,            0,
-           0, cos(alpha), -sin(alpha) ,
-           0, sin(alpha),  cos(alpha);
-      V1 = (R*V1.transpose()).transpose();
-
-      assert(V1(0,2) < 10e-10);
-      assert(V1(1,2) < 10e-10);
-      assert(V1(2,2) < 10e-10);
-
-      // measure the angle between the reference frames
-      // k_ij is the angle between the triangle on the left and the one on the right
-      VectorXd ref0 = V0.row(1) - V0.row(0);
-      VectorXd ref1 = V1.row(1) - V1.row(0);
-
-      ref0.normalize();
-      ref1.normalize();
-
-      double ktemp = atan2(ref1(1),ref1(0)) - atan2(ref0(1),ref0(0));
-
-      // just to be sure, rotate ref0 using angle ktemp...
-      MatrixXd R2(2,2);
-      R2 << cos(ktemp), -sin(ktemp), sin(ktemp), cos(ktemp);
-
-      tmp = R2*ref0.head<2>();
-
-      assert(tmp(0) - ref1(0) < (0.000001));
-      assert(tmp(1) - ref1(1) < (0.000001));
-
-      k[eid] = ktemp;
-    }
-  }
-
-}
-
-
-  void FrameInterpolator::frame2canonical(const Eigen::MatrixXd& TP, const Eigen::RowVectorXd& v, double& theta, Eigen::VectorXd& S_v)
-{
-  using namespace std;
-  using namespace Eigen;
-
-  RowVectorXd v0 = v.segment<3>(0);
-  RowVectorXd v1 = v.segment<3>(3);
-
-  // Project onto the tangent plane
-  Vector2d vp0 = TP * v0.transpose();
-  Vector2d vp1 = TP * v1.transpose();
-
-  // Assemble matrix
-  MatrixXd M(2,2);
-  M << vp0, vp1;
-
-  if (M.determinant() < 0)
-    M.col(1) = -M.col(1);
-
-  assert(M.determinant() > 0);
-
-  // cerr << "M: " << M << endl;
-
-  MatrixXd R,S;
-  PolarDecomposition(M,R,S);
-
-  // Finally, express the cross field as an angle
-  theta = atan2(R(1,0),R(0,0));
-
-  MatrixXd R2(2,2);
-  R2 << cos(theta), -sin(theta), sin(theta), cos(theta);
-
-  assert((R2-R).norm() < 10e-8);
-
-  // Convert into rotation invariant form
-  S = R * S * R.inverse();
-
-  // Copy in vector form
-  S_v = VectorXd(3);
-  S_v << S(0,0), S(0,1), S(1,1);
-}
-
-  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_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 << cos(theta), -sin(theta), sin(theta), cos(theta);
-
-  // Rotation invariant reconstruction
-  MatrixXd M = S_temp * R;
-
-  Vector2d vp0(M(0,0),M(1,0));
-  Vector2d vp1(M(0,1),M(1,1));
-
-  // Unproject the vectors
-  RowVectorXd v0 = vp0.transpose() * TP;
-  RowVectorXd v1 = vp1.transpose() * TP;
-
-  v.resize(6);
-  v << v0, v1;
-}
-
-void FrameInterpolator::solve()
-{
-  interpolateCross();
-  interpolateSymmetric();
-}
-
-void FrameInterpolator::interpolateSymmetric()
-{
-  using namespace std;
-  using namespace Eigen;
-
-  // Generate uniform Laplacian matrix
-  typedef Eigen::Triplet<double> triplet;
-  std::vector<triplet> triplets;
-
-  // Variables are stacked as x1,y1,z1,x2,y2,z2
-  triplets.reserve(3*4*F.rows());
-
-  MatrixXd b = MatrixXd::Zero(3*F.rows(),1);
-
-  // Build L and b
-  for (unsigned eid=0; eid<EF.rows(); ++eid)
-  {
-    if (!isBorderEdge[eid])
-    {
-      for (int z=0;z<2;++z)
-      {
-        // W = [w_a, w_b
-        //      w_b, w_c]
-        //
-
-        // It is not symmetric
-        int i    = EF(eid,z==0?0:1);
-        int j    = EF(eid,z==0?1:0);
-
-        int w_a_0 = (i*3)+0;
-        int w_b_0 = (i*3)+1;
-        int w_c_0 = (i*3)+2;
-
-        int w_a_1 = (j*3)+0;
-        int w_b_1 = (j*3)+1;
-        int w_c_1 = (j*3)+2;
-
-        // Rotation to change frame
-        double r_a =  cos(z==1?k(eid):-k(eid));
-        double r_b = -sin(z==1?k(eid):-k(eid));
-        double r_c =  sin(z==1?k(eid):-k(eid));
-        double r_d =  cos(z==1?k(eid):-k(eid));
-
-        // First term
-        // w_a_0 = r_a^2 w_a_1 + 2 r_a r_b w_b_1 + r_b^2 w_c_1 = 0
-        triplets.push_back(triplet(w_a_0,w_a_0,                -1 ));
-        triplets.push_back(triplet(w_a_0,w_a_1,           r_a*r_a ));
-        triplets.push_back(triplet(w_a_0,w_b_1,       2 * r_a*r_b ));
-        triplets.push_back(triplet(w_a_0,w_c_1,           r_b*r_b ));
-
-        // Second term
-        // w_b_0 = r_a r_c w_a + (r_b r_c + r_a r_d) w_b + r_b r_d w_c
-        triplets.push_back(triplet(w_b_0,w_b_0,                -1 ));
-        triplets.push_back(triplet(w_b_0,w_a_1,           r_a*r_c ));
-        triplets.push_back(triplet(w_b_0,w_b_1, r_b*r_c + r_a*r_d ));
-        triplets.push_back(triplet(w_b_0,w_c_1,           r_b*r_d ));
-
-        // Third term
-        // w_c_0 = r_c^2 w_a + 2 r_c r_d w_b +  r_d^2 w_c
-        triplets.push_back(triplet(w_c_0,w_c_0,                -1 ));
-        triplets.push_back(triplet(w_c_0,w_a_1,           r_c*r_c ));
-        triplets.push_back(triplet(w_c_0,w_b_1,       2 * r_c*r_d ));
-        triplets.push_back(triplet(w_c_0,w_c_1,           r_d*r_d ));
-      }
-    }
-  }
-
-  SparseMatrix<double> L(3*F.rows(),3*F.rows());
-  L.setFromTriplets(triplets.begin(), triplets.end());
-
-  triplets.clear();
-
-  // Add soft constraints
-  double w = 100000;
-  for (unsigned fid=0; fid < F.rows(); ++fid)
-  {
-    if (S_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(fid,i);
-      }
-    }
-  }
-
-  SparseMatrix<double> soft(3*F.rows(),3*F.rows());
-  soft.setFromTriplets(triplets.begin(), triplets.end());
-
-  SparseMatrix<double> M;
-
-  M = L + soft;
-
-  // Solve Lx = b;
-
-  SparseLU<SparseMatrix<double> > solver;
-
-  solver.compute(M);
-
-  if(solver.info()!=Success)
-  {
-    std::cerr << "LU failed - frame_interpolator.cpp" << std::endl;
-    assert(0);
-  }
-
-  MatrixXd x;
-  x = solver.solve(b);
-
-  if(solver.info()!=Success)
-  {
-    std::cerr << "Linear solve failed - frame_interpolator.cpp" << std::endl;
-    assert(0);
-  }
-
-  S = MatrixXd::Zero(F.rows(),3);
-
-  // Copy back the result
-  for (unsigned i=0;i<F.rows();++i)
-    S.row(i) << x(i*3+0), x(i*3+1), x(i*3+2);
-
-}
-
-void FrameInterpolator::setConstraint(const int fid, const Eigen::VectorXd& v)
-{
-  using namespace std;
-  using namespace Eigen;
-
-  double   t_;
-  VectorXd S_;
-
-  frame2canonical(TPs[fid],v,t_,S_);
-
-  Eigen::RowVectorXd v2;
-  canonical2frame(TPs[fid], t_, S_, v2);
-
-  thetas(fid)   = t_;
-  thetas_c[fid] = true;
-
-  S.row(fid) = S_;
-  S_c[fid]   = true;
-
-}
-
-Eigen::MatrixXd FrameInterpolator::getFieldPerFace()
-{
-  using namespace std;
-  using namespace Eigen;
-
-  MatrixXd R(F.rows(),6);
-  for (unsigned i=0; i<F.rows(); ++i)
-  {
-    RowVectorXd v;
-    canonical2frame(TPs[i],thetas(i),S.row(i),v);
-    R.row(i) = v;
-  }
-  return R;
-}
-
-  void FrameInterpolator::PolarDecomposition(Eigen::MatrixXd V, Eigen::MatrixXd& U, Eigen::MatrixXd& P)
-{
-  using namespace std;
-  using namespace Eigen;
-
-  // Polar Decomposition
-  JacobiSVD<MatrixXd> svd(V,Eigen::ComputeFullU | Eigen::ComputeFullV);
-
-  U = svd.matrixU() * svd.matrixV().transpose();
-  P = svd.matrixV() * svd.singularValues().asDiagonal() * svd.matrixV().transpose();
-}
-
-}
-}
-}
-
-IGL_INLINE void igl::copyleft::comiso::frame_field(
-                                 const Eigen::MatrixXd& V,
-                                 const Eigen::MatrixXi& F,
-                                 const Eigen::VectorXi& b,
-                                 const Eigen::MatrixXd& bc1,
-                                 const Eigen::MatrixXd& bc2,
-                                 Eigen::MatrixXd& FF1,
-                                 Eigen::MatrixXd& FF2
-                                 )
-
-{
-  using namespace std;
-  using namespace Eigen;
-
-  assert(b.size() > 0);
-
-  // Init Solver
-  FrameInterpolator field(V,F);
-
-  for (unsigned i=0; i<b.size(); ++i)
-  {
-    VectorXd t(6); t << bc1.row(i).transpose(), bc2.row(i).transpose();
-    field.setConstraint(b(i), t);
-  }
-
-  // Solve
-  field.solve();
-
-  // Copy back
-  MatrixXd R = field.getFieldPerFace();
-  FF1 = R.block(0, 0, R.rows(), 3);
-  FF2 = R.block(0, 3, R.rows(), 3);
-}

include/igl/copyleft/comiso/frame_field.h

@@ -1,54 +0,0 @@
-// This file is part of libigl, a simple c++ geometry processing library.
-//
-// Copyright (C) 2014 Daniele Panozzo <[email protected]>
-//
-// This Source Code Form is subject to the terms of the Mozilla Public License
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can
-// obtain one at http://mozilla.org/MPL/2.0/.
-#ifndef IGL_COMISO_FRAMEFIELD_H
-#define IGL_COMISO_FRAMEFIELD_H
-
-#include "../../igl_inline.h"
-#include "../../PI.h"
-#include <Eigen/Dense>
-#include <vector>
-
-namespace igl
-{
-namespace copyleft
-{
-namespace comiso
-{
-/// Generate a piecewise-constant frame-field field from a sparse set of constraints on faces
-/// using the algorithm proposed in:
-/// Frame Fields: Anisotropic and Non-Orthogonal Cross Fields
-/// Daniele Panozzo, Enrico Puppo, Marco Tarini, Olga Sorkine-Hornung,
-/// ACM Transactions on Graphics (SIGGRAPH, 2014)
-///
-/// @param[in] V       #V by 3 list of mesh vertex coordinates
-/// @param[in] F       #F by 3 list of mesh faces (must be triangles)
-/// @param[in] b       #B by 1 list of constrained face indices
-/// @param[in] bc1     #B by 3 list of the constrained first representative vector of the frame field (up to permutation and sign)
-/// @param[in] bc2     #B by 3 list of the constrained second representative vector of the frame field (up to permutation and sign)
-/// @param[out] FF1      #F by 3 the first representative vector of the frame field (up to permutation and sign)
-/// @param[out] FF2      #F by 3 the second representative vector of the frame field (up to permutation and sign)
-///
-/// \note it now supports only soft constraints, should be extended to support both hard and soft constraints
-IGL_INLINE void frame_field(
-  const Eigen::MatrixXd& V,
-  const Eigen::MatrixXi& F,
-  const Eigen::VectorXi& b,
-  const Eigen::MatrixXd& bc1,
-  const Eigen::MatrixXd& bc2,
-  Eigen::MatrixXd& FF1,
-  Eigen::MatrixXd& FF2
-  );
-}
-}
-}
-
-#ifndef IGL_STATIC_LIBRARY
-#  include "frame_field.cpp"
-#endif
-
-#endif

include/igl/copyleft/comiso/miq.cpp

@@ -1,1543 +0,0 @@
-// This file is part of libigl, a simple c++ geometry processing library.
-//
-// Copyright (C) 2014 Daniele Panozzo <[email protected]>, Olga Diamanti <[email protected]>, Kevin Walliman <[email protected]>
-//
-// 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 "miq.h"
-#include "../../local_basis.h"
-#include "../../triangle_triangle_adjacency.h"
-#include "../../cut_mesh.h"
-#include "../../LinSpaced.h"
-
-// includes for VertexIndexing
-#include "../../HalfEdgeIterator.h"
-#include "../../is_border_vertex.h"
-#include "../../vertex_triangle_adjacency.h"
-
-// includes for PoissonSolver
-#include "../../slice_into.h"
-#include "../../grad.h"
-#include "../../cotmatrix.h"
-#include "../../doublearea.h"
-#include <gmm/gmm.h>
-#include <CoMISo/Solver/ConstrainedSolver.hh>
-#include <CoMISo/Solver/MISolver.hh>
-#include <CoMISo/Solver/GMM_Tools.hh>
-
-//
-#include "igl/cross_field_mismatch.h"
-#include "../../comb_frame_field.h"
-#include "../../comb_cross_field.h"
-#include "../../cut_mesh_from_singularities.h"
-#include "../../find_cross_field_singularities.h"
-#include "../../compute_frame_field_bisectors.h"
-#include "../../rotate_vectors.h"
-
-#ifndef NDEBUG
-#include <fstream>
-#endif
-#include <iostream>
-#include "../../matlab_format.h"
-
-#define DEBUGPRINT 0
-
-
-namespace igl {
-namespace copyleft {
-namespace comiso {
-  struct SeamInfo
-  {
-    int v0,v0p;
-    unsigned int integerVar;
-    int mismatch;
-
-    IGL_INLINE SeamInfo(int _v0,
-                        int _v0p,
-                        int _mismatch,
-                        unsigned int _integerVar);
-
-    IGL_INLINE SeamInfo(const SeamInfo &S1);
-  };
-
-  struct MeshSystemInfo
-  {
-    MeshSystemInfo()
-    {
-      num_vert_variables = 0;
-      num_integer_cuts = 0;
-    }
-    ////number of vertices variables
-    unsigned int num_vert_variables;
-    ///num of integer for cuts
-    unsigned int num_integer_cuts;
-    ///this are used for drawing purposes
-    std::vector<SeamInfo> edgeSeamInfo;
-  };
-
-
-  template <typename DerivedV, typename DerivedF>
-  class VertexIndexing
-  {
-  public:
-    // Input:
-    const Eigen::PlainObjectBase<DerivedV> &V;
-    const Eigen::PlainObjectBase<DerivedF> &F;
-    const Eigen::PlainObjectBase<DerivedV> &Vcut;
-    const Eigen::PlainObjectBase<DerivedF> &Fcut;
-    const Eigen::PlainObjectBase<DerivedF> &TT;
-    const Eigen::PlainObjectBase<DerivedF> &TTi;
-
-    const Eigen::Matrix<int, Eigen::Dynamic, 3> &mismatch;
-    const Eigen::Matrix<int, Eigen::Dynamic, 1> &singular; // bool
-    const Eigen::Matrix<int, Eigen::Dynamic, 3> &seams; // 3 bool
-
-
-    ///this handle for mesh TODO: move with the other global variables
-    MeshSystemInfo systemInfo;
-
-    IGL_INLINE VertexIndexing(const Eigen::PlainObjectBase<DerivedV> &_V,
-                              const Eigen::PlainObjectBase<DerivedF> &_F,
-                              const Eigen::PlainObjectBase<DerivedV> &_Vcut,
-                              const Eigen::PlainObjectBase<DerivedF> &_Fcut,
-                              const Eigen::PlainObjectBase<DerivedF> &_TT,
-                              const Eigen::PlainObjectBase<DerivedF> &_TTi,
-                              const Eigen::Matrix<int, Eigen::Dynamic, 3> &_mismatch,
-                              const Eigen::Matrix<int, Eigen::Dynamic, 1> &_singular,
-                              const Eigen::Matrix<int, Eigen::Dynamic, 3> &_seams
-                              );
-
-    // provide information about every vertex per seam
-    IGL_INLINE void initSeamInfo();
-
-
-  private:
-    struct VertexInfo{
-      int v;  // vertex index (according to V)
-      int f0, k0; // face and local edge information of the edge that connects this vertex to the previous vertex (previous in the vector)
-      int f1, k1; // face and local edge information of the other face corresponding to the same edge
-      VertexInfo(int _v, int _f0, int _k0, int _f1, int _k1) :
-                 v(_v), f0(_f0), k0(_k0), f1(_f1), k1(_k1){}
-      bool operator==(VertexInfo const& other){
-        return other.v == v;
-      }
-    };
-
-    IGL_INLINE void getSeamInfo(int f0,
-                                int f1,
-                                int indexE,
-                                int &v0, int &v1,
-                                int &v0p, int &v1p,
-                                int &_mismatch);
-
-    IGL_INLINE std::vector<std::vector<VertexInfo> > getVerticesPerSeam();
-  };
-
-
-  template <typename DerivedV, typename DerivedF>
-  class PoissonSolver
-  {
-  private:
-
-    // Penalization term for integer variables used in mixedIntegerSolve
-    const double PENALIZATION = 0.000001;
-
-  public:
-    IGL_INLINE void solvePoisson(Eigen::VectorXd stiffness,
-                                 double gradientSize = 0.1,
-                                 double gridResolution = 1.,
-                                 bool directRound = true,
-                                 unsigned int localIter = 0,
-                                 bool doRound = true,
-                                 bool singularityRound = true,
-                                 const std::vector<int> &roundVertices = std::vector<int>(),
-                                 const 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<DerivedV> &_Vcut,
-                             const Eigen::PlainObjectBase<DerivedF> &_Fcut,
-                             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>&_singular,
-                             const MeshSystemInfo &_systemInfo
-                             );
-
-    const Eigen::PlainObjectBase<DerivedV> &V;
-    const Eigen::PlainObjectBase<DerivedF> &F;
-    const Eigen::PlainObjectBase<DerivedV> &Vcut;
-    const Eigen::PlainObjectBase<DerivedF> &Fcut;
-    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> &singular; // bool
-
-    const MeshSystemInfo &systemInfo;
-
-    // Internal:
-    Eigen::VectorXd Handle_Stiffness;
-    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;
-    // per vertex UV coordinates, Vcut.rows() x 2
-    Eigen::MatrixXd UV_out;
-
-    // Matrices
-    Eigen::SparseMatrix<double> Lhs;
-    Eigen::SparseMatrix<double> Constraints;
-    Eigen::VectorXd rhs;
-    Eigen::VectorXd constraints_rhs;
-    ///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;
-
-    ///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 COMMON MATH FUNCTIONS
-    ///return the complex encoding the rotation
-    ///for a given mismatch interval
-    IGL_INLINE std::complex<double> getRotationComplex(int interval);
-    ///END COMMON MATH FUNCTIONS
-
-    ///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 cycling over all range maps
-    //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 coneGridRes = 1,
-                                      bool directRound = true,
-                                      unsigned 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;
-    DerivedV Vcut;
-    DerivedF Fcut;
-    Eigen::MatrixXd UV_out;
-    DerivedF FUV_out;
-
-    // internal
-    DerivedF TT;
-    DerivedF TTi;
-
-    // Stiffness per face
-    Eigen::VectorXd stiffnessVector;
-    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::Matrix<int, Eigen::Dynamic, 3> &mismatch,
-                         const Eigen::Matrix<int, Eigen::Dynamic, 1> &singular,
-                         const Eigen::Matrix<int, Eigen::Dynamic, 3> &seams,
-                         Eigen::PlainObjectBase<DerivedU> &UV,
-                         Eigen::PlainObjectBase<DerivedF> &FUV,
-                         double gradientSize = 30.0,
-                         double stiffness = 5.0,
-                         bool directRound = false,
-                         unsigned int iter = 5,
-                         unsigned int localIter = 5,
-                         bool doRound = true,
-                         bool singularityRound = 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(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(int i, const Eigen::MatrixXd& WUV);
-
-  };
-};
-};
-}
-
-IGL_INLINE igl::copyleft::comiso::SeamInfo::SeamInfo(int _v0,
-                                                     int _v0p,
-                                                     int _mismatch,
-                                                     unsigned int _integerVar)
-{
-  v0=_v0;
-  v0p=_v0p;
-  integerVar=_integerVar;
-  mismatch=_mismatch;
-}
-
-IGL_INLINE igl::copyleft::comiso::SeamInfo::SeamInfo(const SeamInfo &S1)
-{
-  v0=S1.v0;
-  v0p=S1.v0p;
-  integerVar=S1.integerVar;
-  mismatch=S1.mismatch;
-}
-
-
-template <typename DerivedV, typename DerivedF>
-IGL_INLINE igl::copyleft::comiso::VertexIndexing<DerivedV, DerivedF>::VertexIndexing(const Eigen::PlainObjectBase<DerivedV> &_V,
-                                                                   const Eigen::PlainObjectBase<DerivedF> &_F,
-                                                                   const Eigen::PlainObjectBase<DerivedV> &_Vcut,
-                                                                   const Eigen::PlainObjectBase<DerivedF> &_Fcut,
-                                                                   const Eigen::PlainObjectBase<DerivedF> &_TT,
-                                                                   const Eigen::PlainObjectBase<DerivedF> &_TTi,
-                                                                   const Eigen::Matrix<int, Eigen::Dynamic, 3> &_mismatch,
-                                                                   const Eigen::Matrix<int, Eigen::Dynamic, 1> &_singular,
-                                                                   const Eigen::Matrix<int, Eigen::Dynamic, 3> &_seams
-                                                                   ):
-V(_V),
-F(_F),
-Vcut(_Vcut),
-Fcut(_Fcut),
-TT(_TT),
-TTi(_TTi),
-mismatch(_mismatch),
-singular(_singular),
-seams(_seams)
-{
-  #ifdef DEBUG_PRINT
-  cerr<<igl::matlab_format(Handle_Seams,"Handle_Seams");
-#endif
-
-  systemInfo.num_vert_variables=Vcut.rows();
-  systemInfo.num_integer_cuts=0;
-}
-
-template <typename DerivedV, typename DerivedF>
-IGL_INLINE void igl::copyleft::comiso::VertexIndexing<DerivedV, DerivedF>::getSeamInfo(const int f0,
-                                                                                       const int f1,
-                                                                                       const int indexE,
-                                                                                       int &v0, int &v1,
-                                                                                       int &v0p, int &v1p,
-                                                                                       int &_mismatch)
-{
-  int edgef0 = indexE;
-  v0 = Fcut(f0,edgef0);
-  v1 = Fcut(f0,(edgef0+1)%3);
-  ////get the index on opposite side
-  assert(TT(f0,edgef0) == f1);
-  int edgef1 = TTi(f0,edgef0);
-  v1p = Fcut(f1,edgef1);
-  v0p = Fcut(f1,(edgef1+1)%3);
-
-  _mismatch = mismatch(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 std::vector<std::vector<typename igl::copyleft::comiso::VertexIndexing<DerivedV, DerivedF>::VertexInfo> > igl::copyleft::comiso::VertexIndexing<DerivedV, DerivedF>::getVerticesPerSeam()
-{
-  // Return value
-  std::vector<std::vector<VertexInfo> >verticesPerSeam;
-
-  // for every vertex, keep track of their adjacent vertices on seams.
-  // regular vertices have two neighbors on a seam, start- and endvertices may have any other numbers of neighbors (e.g. 1 or 3)
-  std::vector<std::list<VertexInfo> > VVSeam(V.rows());
-  Eigen::MatrixXi F_hit = Eigen::MatrixXi::Zero(F.rows(), 3);
-  for (unsigned int f=0; f<F.rows();f++)
-  {
-    int f0 = f;
-    for(int k0=0; k0<3; k0++){
-      int f1 = TT(f0,k0);
-      if(f1 == -1)
-        continue;
-
-      if (seams(f0, k0) != 0 && F_hit(f0,k0) == 0)
-      {
-        int v0 = F(f0, k0);
-        int v1 = F(f0, (k0+1)%3);
-        int k1 = TTi(f0,k0);
-        VVSeam[v0].push_back(VertexInfo(v1, f0, k0, f1, k1));
-        VVSeam[v1].push_back(VertexInfo(v0, f0, k0, f1, k1));
-        F_hit(f0, k0) = 1;
-        F_hit(f1, k1) = 1;
-      }
-    }
-  }
-
-  // Find start vertices, i.e. vertices that start or end a seam branch
-  std::vector<int> startVertexIndices;
-  std::vector<bool> isStartVertex(V.rows());
-  for (unsigned int i=0;i<V.rows();i++)
-  {
-    isStartVertex[i] = false;
-    // vertices with two neighbors are regular vertices, unless the vertex is a singularity, in which case it qualifies as a start vertex
-    if ((!VVSeam[i].empty() && VVSeam[i].size() != 2) || singular(i) != 0)
-    {
-      startVertexIndices.push_back(i);
-      isStartVertex[i] = true;
-    }
-  }
-
-  // For each startVertex, walk along its seam
-  for (auto element : startVertexIndices)
-  {
-    auto startVertexNeighbors = &VVSeam[element];
-    size_t neighborSize = startVertexNeighbors->size();
-
-    // explore every seam to which this vertex is a start vertex
-    // note: a vertex can never be a start vertex and a regular vertex simultaneously
-    for (size_t j = 0; j < neighborSize; j++)
-    {
-      std::vector<VertexInfo> thisSeam; // temporary container
-
-      // Create vertexInfo struct for start vertex
-      VertexInfo startVertex = VertexInfo(element, -1, -1, -1, -1);// -1 values are arbitrary (will never be used)
-      VertexInfo currentVertex = startVertex;
-      // Add start vertex to the seam
-      thisSeam.push_back(currentVertex);
-
-      // advance on the seam
-      auto currentVertexNeighbors = startVertexNeighbors;
-      auto nextVertex = currentVertexNeighbors->front();
-      currentVertexNeighbors->pop_front();
-
-      // bogus initialization due to lack of def. constructor
-      VertexInfo prevVertex = startVertex;
-      while (true)
-      {
-        // move to the next vertex
-        prevVertex = currentVertex;
-        currentVertex = nextVertex;
-        currentVertexNeighbors = &VVSeam[nextVertex.v];
-
-        // add current vertex to this seam
-        thisSeam.push_back(currentVertex);
-
-        // remove the previous vertex
-        auto it = std::find(currentVertexNeighbors->begin(), currentVertexNeighbors->end(), prevVertex);
-        assert(it != currentVertexNeighbors->end());
-        currentVertexNeighbors->erase(it);
-
-        if (currentVertexNeighbors->size() == 1 && !isStartVertex[currentVertex.v])
-        {
-          nextVertex = currentVertexNeighbors->front();
-          currentVertexNeighbors->pop_front();
-        }
-        else
-          break;
-      }
-      verticesPerSeam.push_back(thisSeam);
-    }
-  }
-
-  return verticesPerSeam;
-}
-
-template <typename DerivedV, typename DerivedF>
-IGL_INLINE void igl::copyleft::comiso::VertexIndexing<DerivedV, DerivedF>::initSeamInfo()
-{
-  auto verticesPerSeam = getVerticesPerSeam();
-  systemInfo.edgeSeamInfo.clear();
-  unsigned int integerVar = 0;
-  // Loop over each seam
-  for(auto seam : verticesPerSeam){
-    //choose initial side of the seam such that the start vertex corresponds to Fcut(f, k) and the end vertex corresponds to Fcut(f, (k+1)%3) and not vice versa.
-    int priorVertexIdx;
-    if(seam.size() > 2){
-      auto v1 = seam[1];
-      auto v2 = seam[2];
-      if(Fcut(v1.f0, (v1.k0+1) % 3) == Fcut(v2.f0, v2.k0) || Fcut(v1.f0, (v1.k0+1) % 3) == Fcut(v2.f1, v2.k1)){
-        priorVertexIdx = Fcut(v1.f0, v1.k0);
-      }
-      else{
-        priorVertexIdx = Fcut(v1.f1, v1.k1);
-        assert(Fcut(v1.f1, (v1.k1+1) % 3) == Fcut(v2.f0, v2.k0) || Fcut(v1.f1, (v1.k1+1) % 3) == Fcut(v2.f1, v2.k1));
-      }
-    }
-    else{
-      auto v1 = seam[1];
-      priorVertexIdx = Fcut(v1.f0, v1.k0);
-    }
-
-    // Loop over each vertex of the seam
-    for(auto it=seam.begin()+1; it != seam.end(); ++it){
-      auto vertex = *it;
-      // choose the correct side of the seam
-      int f,k,ff;
-      if(priorVertexIdx == Fcut(vertex.f0, vertex.k0)){
-        f = vertex.f0; ff = vertex.f1;
-        k = vertex.k0;
-      }
-      else{
-        f = vertex.f1; ff = vertex.f0;
-        k = vertex.k1;
-        assert(priorVertexIdx == Fcut(vertex.f1, vertex.k1));
-      }
-
-      int vtx0,vtx0p,vtx1,vtx1p;
-      int MM;
-      getSeamInfo(f, ff, k, vtx0, vtx1, vtx0p, vtx1p, MM);
-      systemInfo.edgeSeamInfo.push_back(SeamInfo(vtx0,vtx0p,MM,integerVar));
-      if(it == seam.end() -1){
-        systemInfo.edgeSeamInfo.push_back(SeamInfo(vtx1,vtx1p,MM,integerVar));
-      }
-      priorVertexIdx = vtx1;
-    }
-    // use the same integer for each seam
-    integerVar++;
-  }
-  systemInfo.num_integer_cuts = integerVar;
-
-#ifndef NDEBUG
-  int totalNVerticesOnSeams = 0;
-  for(auto const & seam : verticesPerSeam){
-    totalNVerticesOnSeams += seam.size();
-  }
-  assert(systemInfo.edgeSeamInfo.size() == totalNVerticesOnSeams);
-#endif
-}
-
-
-
-template <typename DerivedV, typename DerivedF>
-IGL_INLINE void igl::copyleft::comiso::PoissonSolver<DerivedV, DerivedF>::solvePoisson(Eigen::VectorXd stiffness,
-                                                                                       double gradientSize,
-                                                                                       double gridResolution,
-                                                                                       bool directRound,
-                                                                                       unsigned int localIter,
-                                                                                       bool doRound,
-                                                                                       bool singularityRound,
-                                                                                       const std::vector<int> &roundVertices,
-                                                                                       const std::vector<std::vector<int>> &hardFeatures)
-{
-  Handle_Stiffness = stiffness;
-
-  //initialization of flags and data structures
-  integer_rounding=doRound;
-
-  ids_to_round.clear();
-
-  clearUserConstraint();
-  // copy the user constraints number
-  for (const auto & element : hardFeatures)
-  {
-    addSharpEdgeConstraint(element[0], element[1]);
-  }
-
-  ///Initializing Matrix
-  clock_t t0 = clock();
-
-  ///initialize the matrix ALLOCATING SPACE
-  initMatrix();
-  if (DEBUGPRINT)
-    printf("\n ALLOCATED THE MATRIX \n");
-
-  ///build the Laplacian system
-  buildLaplacianMatrix(gradientSize);
-
-  // 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)
-    addToRoundVertices(roundVertices);
-
-  if (singularityRound)
-    addSingularityRound();
-
-  clock_t t1 = clock();
-  if (DEBUGPRINT) printf("\n time:%ld \n",t1-t0);
-  if (DEBUGPRINT) printf("\n SOLVING \n");
-
-  mixedIntegerSolve(gridResolution, directRound, localIter);
-
-  clock_t t2 = clock();
-  if (DEBUGPRINT) printf("\n time:%ld \n",t2-t1);
-  if (DEBUGPRINT) printf("\n ASSIGNING COORDS \n");
-
-  mapCoords();
-
-  clock_t t3 = clock();
-  if (DEBUGPRINT) printf("\n time:%ld \n",t3-t2);
-  if (DEBUGPRINT) printf("\n FINISHED \n");
-}
-
-template <typename DerivedV, typename DerivedF>
-IGL_INLINE igl::copyleft::comiso::PoissonSolver<DerivedV, DerivedF>
-::PoissonSolver(const Eigen::PlainObjectBase<DerivedV> &_V,
-                const Eigen::PlainObjectBase<DerivedF> &_F,
-                const Eigen::PlainObjectBase<DerivedV> &_Vcut,
-                const Eigen::PlainObjectBase<DerivedF> &_Fcut,
-                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>&_singular,
-                const MeshSystemInfo &_systemInfo
-):
-V(_V),
-F(_F),
-Vcut(_Vcut),
-Fcut(_Fcut),
-TT(_TT),
-TTi(_TTi),
-PD1(_PD1),
-PD2(_PD2),
-singular(_singular),
-systemInfo(_systemInfo)
-{
-  n_fixed_vars = 0;
-  n_vert_vars = 0;
-  num_total_vars = 0;
-  n_integer_vars = 0;
-  num_cut_constraint = 0;
-  num_userdefined_constraint = 0;
-  num_constraint_equations = 0;
-  integer_rounding = false;
-  UV = Eigen::MatrixXd(V.rows(),2);
-  WUV = Eigen::MatrixXd(F.rows(),6);
-  UV_out = Eigen::MatrixXd(Vcut.rows(),2);
-  igl::vertex_triangle_adjacency(V,F,VF,VFi);
-}
-
-///START COMMON MATH FUNCTIONS
-///return the complex encoding the rotation
-///for a given mismatch interval
-template <typename DerivedV, typename DerivedF>
-IGL_INLINE std::complex<double> igl::copyleft::comiso::PoissonSolver<DerivedV, DerivedF>::getRotationComplex(
-        int interval)
-{
-  assert((interval>=0)&&(interval<4));
-
-  switch(interval)
-  {
-    case 0:return {1,0};
-    case 1:return {0,1};
-    case 2: return {-1,0};
-    default:return {0,-1};
-  }
-}
-
-///END COMMON MATH FUNCTIONS
-
-///START FIXING VERTICES
-///set a given vertex as fixed
-template <typename DerivedV, typename DerivedF>
-IGL_INLINE void igl::copyleft::comiso::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::copyleft::comiso::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 (singular(v))
-    {
-      addFixedVertex(v);
-      UV.row(v) << 0,0;
-      return;
-    }
-  }
-
-  ///if anything fixed fix the first
-  addFixedVertex(0);
-  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::copyleft::comiso::PoissonSolver<DerivedV, DerivedF>::findFixedVert()
-{
-  Hard_constraints.clear();
-  findFixedVertField();
-}
-
-template <typename DerivedV, typename DerivedF>
-IGL_INLINE int igl::copyleft::comiso::PoissonSolver<DerivedV, DerivedF>::getFirstVertexIndex(int v)
-{
-  return Fcut(VF[v][0],VFi[v][0]);
-}
-
-///fix the vertices which are flagged as fixed
-template <typename DerivedV, typename DerivedF>
-IGL_INLINE void igl::copyleft::comiso::PoissonSolver<DerivedV, DerivedF>::fixBlockedVertex()
-{
-  int offset_row = 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
-    Constraints.coeffRef(indexRow,  indexvert)   += 1;
-    Constraints.coeffRef(indexRow+1,indexvert+1) += 1;
-
-    ///add fixing constraint RHS
-    constraints_rhs[indexCol]   = UV(v,0);
-    constraints_rhs[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::copyleft::comiso::PoissonSolver<DerivedV, DerivedF>::addSingularityRound()
-{
-  for (unsigned int v=0;v<V.rows();v++)
-  {
-    if (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::copyleft::comiso::PoissonSolver<DerivedV, DerivedF>::addToRoundVertices(std::vector<int> ids)
-{
-  for(auto index : ids)
-  {
-    if (index < 0 || index >= V.rows())
-      std::cerr << "WARNING: Ignored round vertex constraint, vertex " << index << " does not exist in the mesh." << std::endl;
-    int index0 = getFirstVertexIndex(index);
-    ids_to_round.push_back( index0 * 2   );
-    ids_to_round.push_back((index0 * 2)+1);
-  }
-}
-
-///START GENERIC SYSTEM FUNCTIONS
-template <typename DerivedV, typename DerivedF>
-IGL_INLINE void igl::copyleft::comiso::PoissonSolver<DerivedV, DerivedF>::buildLaplacianMatrix(double vfscale)
-{
-  Eigen::VectorXi idx  = igl::LinSpaced<Eigen::VectorXi >(Vcut.rows(), 0, 2*Vcut.rows()-2);
-  Eigen::VectorXi idx2 = igl::LinSpaced<Eigen::VectorXi >(Vcut.rows(), 1, 2*Vcut.rows()-1);
-
-  // get gradient matrix
-  Eigen::SparseMatrix<double> G(Fcut.rows() * 3, Vcut.rows());
-  igl::grad(Vcut, Fcut, G);
-
-  // get triangle weights
-  Eigen::VectorXd dblA(Fcut.rows());
-  igl::doublearea(Vcut, Fcut, dblA);
-
-  // compute intermediate result
-  Eigen::SparseMatrix<double> G2;
-  G2 = G.transpose() * dblA.replicate<3,1>().asDiagonal() * Handle_Stiffness.replicate<3,1>().asDiagonal();
-
-  ///  Compute LHS
-  Eigen::SparseMatrix<double> Cotmatrix;
-  Cotmatrix = 0.5 * G2 * G;
-  igl::slice_into(Cotmatrix, idx,  idx,  Lhs);
-  igl::slice_into(Cotmatrix, idx2, idx2, Lhs);
-
-  /// Compute RHS
-  // reshape nrosy vectors
-  const Eigen::MatrixXd u = Eigen::Map<const Eigen::MatrixXd>(PD1.data(),Fcut.rows()*3,1); // this mimics a reshape at the cost of a copy.
-  const Eigen::MatrixXd v = Eigen::Map<const Eigen::MatrixXd>(PD2.data(),Fcut.rows()*3,1); // this mimics a reshape at the cost of a copy.
-
-  // multiply with weights
-  Eigen::VectorXd rhs1 =  G2 * u * 0.5 * vfscale;
-  Eigen::VectorXd rhs2 = -G2 * v * 0.5 * vfscale;
-  rhs( idx) = rhs1;
-  rhs(idx2) = rhs2;
-}
-
-///find different sized of the system
-template <typename DerivedV, typename DerivedF>
-IGL_INLINE void igl::copyleft::comiso::PoissonSolver<DerivedV, DerivedF>::findSizes()
-{
-  ///find the vertex that need to be fixed
-  findFixedVert();
-
-  ///REAL PART
-  n_vert_vars = systemInfo.num_vert_variables;
-
-  ///INTEGER PART
-  ///the total number of integer variables
-  n_integer_vars = systemInfo.num_integer_cuts;
-
-  ///CONSTRAINT PART
-  num_cut_constraint = systemInfo.edgeSeamInfo.size();
-
-  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+n_integer_vars) * 2;
-
-  ///initialize matrix size
-
-  if (DEBUGPRINT)     printf("\n*** SYSTEM VARIABLES *** \n");
-  if (DEBUGPRINT)     printf("* NUM REAL VERTEX VARIABLES %ud \n",n_vert_vars);
-
-  if (DEBUGPRINT)     printf("\n*** INTEGER VARIABLES *** \n");
-  if (DEBUGPRINT)     printf("* NUM INTEGER VARIABLES %ud \n",n_integer_vars);
-
-  if (DEBUGPRINT)     printf("\n*** CONSTRAINTS *** \n ");
-  if (DEBUGPRINT)     printf("* NUM FIXED CONSTRAINTS %ud\n",n_fixed_vars);
-  if (DEBUGPRINT)     printf("* NUM CUTS CONSTRAINTS %ud\n",num_cut_constraint);
-  if (DEBUGPRINT)     printf("* NUM USER DEFINED CONSTRAINTS %ud\n",num_userdefined_constraint);
-
-  if (DEBUGPRINT)     printf("\n*** TOTAL SIZE *** \n");
-  if (DEBUGPRINT)     printf("* TOTAL VARIABLE SIZE (WITH INTEGER TRASL) %ud \n",num_total_vars);
-  if (DEBUGPRINT)     printf("* TOTAL CONSTRAINTS %ud \n",num_constraint_equations);
-}
-
-template <typename DerivedV, typename DerivedF>
-IGL_INLINE void igl::copyleft::comiso::PoissonSolver<DerivedV, DerivedF>::allocateSystem()
-{
-  Lhs.resize(n_vert_vars * 2, n_vert_vars * 2);
-  Constraints.resize(num_constraint_equations, num_total_vars);
-  rhs.resize(n_vert_vars * 2);
-  constraints_rhs.resize(num_constraint_equations);
-
-  printf("\n INITIALIZED SPARSE MATRIX OF %ud x %ud \n",n_vert_vars*2, n_vert_vars*2);
-  printf("\n INITIALIZED SPARSE MATRIX OF %ud x %ud \n",num_constraint_equations, num_total_vars);
-  printf("\n INITIALIZED VECTOR OF %ud x 1 \n",n_vert_vars*2);
-  printf("\n INITIALIZED VECTOR OF %ud x 1 \n",num_constraint_equations);
-}
-
-///intitialize the whole matrix
-template <typename DerivedV, typename DerivedF>
-IGL_INLINE void igl::copyleft::comiso::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::copyleft::comiso::PoissonSolver<DerivedV, DerivedF>::mapCoords()
-{
-  ///map coords to faces
-  for (unsigned int f=0;f<Fcut.rows();f++)
-  {
-
-    for (int k=0;k<3;k++)
-    {
-      //get the index of the variable in the system
-      int indexUV = Fcut(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;
-    }
-
-  }
-
-  for(unsigned int i = 0; i < Vcut.rows(); i++){
-    UV_out(i,0) = X[i*2];
-    UV_out(i,1) = X[i*2+1];
-  }
-}
-
-///END GENERIC SYSTEM FUNCTIONS
-
-///set the constraints for the inter-range cuts
-template <typename DerivedV, typename DerivedF>
-IGL_INLINE void igl::copyleft::comiso::PoissonSolver<DerivedV, DerivedF>::buildSeamConstraintsExplicitTranslation()
-{
-  ///current constraint row
-  int constr_row = 0;
-
-  for (unsigned int i=0; i<num_cut_constraint; i++)
-  {
-    int interval = systemInfo.edgeSeamInfo[i].mismatch;
-    if (interval==1)
-      interval=3;
-    else
-      if(interval==3)
-        interval=1;
-
-    int p0  = systemInfo.edgeSeamInfo[i].v0;
-    int p0p = systemInfo.edgeSeamInfo[i].v0p;
-
-    std::complex<double> rot = getRotationComplex(interval);
-
-    ///get the integer variable
-    unsigned int integerVar = n_vert_vars + systemInfo.edgeSeamInfo[i].integerVar;
-
-    if (integer_rounding)
-    {
-      ids_to_round.push_back(integerVar*2);
-      ids_to_round.push_back(integerVar*2+1);
-    }
-
-    // cross boundary compatibility conditions
-    Constraints.coeffRef(constr_row,   2*p0)   +=  rot.real();
-    Constraints.coeffRef(constr_row,   2*p0+1) += -rot.imag();
-    Constraints.coeffRef(constr_row+1, 2*p0)   +=  rot.imag();
-    Constraints.coeffRef(constr_row+1, 2*p0+1) +=  rot.real();
-
-    Constraints.coeffRef(constr_row,   2*p0p)   += -1;
-    Constraints.coeffRef(constr_row+1, 2*p0p+1) += -1;
-
-    Constraints.coeffRef(constr_row,   2*integerVar)   += 1;
-    Constraints.coeffRef(constr_row+1, 2*integerVar+1) += 1;
-
-    constraints_rhs[constr_row]   = 0;
-    constraints_rhs[constr_row+1] = 0;
-
-    constr_row += 2;
-  }
-
-}
-
-///set the constraints for the inter-range cuts
-template <typename DerivedV, typename DerivedF>
-IGL_INLINE void igl::copyleft::comiso::PoissonSolver<DerivedV, DerivedF>::buildUserDefinedConstraints()
-{
-  /// the user defined constraints are at the end
-  unsigned int constr_row = num_cut_constraint*2 + n_fixed_vars*2;
-
-  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)
-    {
-      Constraints.coeffRef(constr_row, j) = userdefined_constraints[i][j];
-    }
-
-    constraints_rhs[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::copyleft::comiso::PoissonSolver<DerivedV, DerivedF>::mixedIntegerSolve(double coneGridRes,
-                                                                                            bool directRound,
-                                                                                            unsigned int localIter)
-{
-  X = std::vector<double>((n_vert_vars+n_integer_vars)*2);
-  if (DEBUGPRINT)
-    printf("\n ALLOCATED X \n");
-
-  ///variables part
-  const int sizeMatrix = (n_vert_vars + n_integer_vars) * 2;
-  const int scalarSize = n_vert_vars * 2;
-
-  ///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> B(sizeMatrix,0);  // rhs
-
-  if (DEBUGPRINT)
-    printf("\n ALLOCATED RHS STRUCTURES \n");
-
-  //// copy LHS
-  for (int k=0; k < Lhs.outerSize(); ++k){
-    for (Eigen::SparseMatrix<double>::InnerIterator it(Lhs,k); it; ++it){
-      int row = it.row();
-      int col = it.col();
-      A(row, col) += it.value();
-    }
-  }
-  //// copy Constraints
-  for (int k=0; k < Constraints.outerSize(); ++k){
-    for (Eigen::SparseMatrix<double>::InnerIterator it(Constraints,k); it; ++it){
-      int row = it.row();
-      int col = it.col();
-      C(row, col) += it.value();
-    }
-  }
-
-  if (DEBUGPRINT)
-    printf("\n SET %d INTEGER VALUES \n",n_integer_vars);
-
-  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(unsigned int i = 0; i < scalarSize; ++i)
-  {
-    B[i] = rhs[i] * coneGridRes;
-  }
-
-  // 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) = -constraints_rhs[i] * coneGridRes;
-  }
-
-  COMISO::ConstrainedSolver solver;
-
-  solver.misolver().set_local_iters(localIter);
-
-  solver.misolver().set_direct_rounding(directRound);
-
-  std::sort(ids_to_round.begin(),ids_to_round.end());
-  auto new_end=std::unique(ids_to_round.begin(),ids_to_round.end());
-  long int dist = distance(ids_to_round.begin(),new_end);
-  ids_to_round.resize(dist);
-
-  solver.solve( C, A, X, B, ids_to_round, 0.0, false, false);
-}
-
-template <typename DerivedV, typename DerivedF>
-IGL_INLINE void igl::copyleft::comiso::PoissonSolver<DerivedV, DerivedF>::clearUserConstraint()
-{
-  num_userdefined_constraint = 0;
-  userdefined_constraints.clear();
-}
-
-template <typename DerivedV, typename DerivedF>
-IGL_INLINE void igl::copyleft::comiso::PoissonSolver<DerivedV, DerivedF>::addSharpEdgeConstraint(int fid, int vid)
-{
-  // prepare constraint
-  std::vector<int> c(systemInfo.num_vert_variables*2 + 1, 0);
-
-  int v1 = Fcut(fid,vid);
-  int v2 = Fcut(fid,(vid+1)%3);
-
-  Eigen::Matrix<typename DerivedV::Scalar, 3, 1> e = Vcut.row(v2) - Vcut.row(v1);
-  e = e.normalized();
-
-  double d1 = fabs(e.dot(PD1.row(fid).normalized()));
-  double d2 = fabs(e.dot(PD2.row(fid).normalized()));
-
-  int offset = 0;
-
-  if (d1>d2)
-    offset = 1;
-
-  ids_to_round.push_back((v1 * 2) + offset);
-  ids_to_round.push_back((v2 * 2) + offset);
-
-  // add constraint
-  c[(v1 * 2) + offset] =  1;
-  c[(v2 * 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::copyleft::comiso::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::Matrix<int, Eigen::Dynamic, 3> &mismatch,
-        const Eigen::Matrix<int, Eigen::Dynamic, 1> &singular,
-        const Eigen::Matrix<int, Eigen::Dynamic, 3> &seams,
-        Eigen::PlainObjectBase<DerivedU> &UV,
-        Eigen::PlainObjectBase<DerivedF> &FUV,
-        double gradientSize,
-        double stiffness,
-        bool directRound,
-        unsigned int iter,
-        unsigned int localIter,
-        bool doRound,
-        bool singularityRound,
-        std::vector<int> roundVertices,
-        std::vector<std::vector<int> > hardFeatures):
-V(V_),
-F(F_)
-{
-  igl::cut_mesh(V, F, seams, Vcut, Fcut);
-
-  igl::local_basis(V,F,B1,B2,B3);
-  igl::triangle_triangle_adjacency(F,TT,TTi);
-
-  // Prepare indexing for the linear system
-  VertexIndexing<DerivedV, DerivedF> VInd(V, F, Vcut, Fcut, TT, TTi, mismatch, singular, seams);
-
-  VInd.initSeamInfo();
-
-  // Assemble the system and solve
-  PoissonSolver<DerivedV, DerivedF> PSolver(V,
-                                            F,
-                                            Vcut,
-                                            Fcut,
-                                            TT,
-                                            TTi,
-                                            PD1_combed,
-                                            PD2_combed,
-                                            singular,
-                                            VInd.systemInfo);
-  stiffnessVector = Eigen::VectorXd::Constant(F.rows(),1);
-
-
-  if (iter > 0) // do stiffening
-  {
-    for (unsigned int i=0;i<iter;i++)
-    {
-      PSolver.solvePoisson(stiffnessVector, gradientSize, 1.f, directRound, localIter, doRound, singularityRound,
-                           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(stiffnessVector, gradientSize, 1.f, directRound, localIter, doRound, singularityRound,
-                         roundVertices, hardFeatures);
-  }
-
-  int nflips=NumFlips(PSolver.WUV);
-  printf("**** END OPTIMIZING #FLIPS %d  ****\n",nflips);
-
-  UV_out = PSolver.UV_out;
-  FUV_out = PSolver.Fcut;
-  fflush(stdout);
-}
-
-template <typename DerivedV, typename DerivedF, typename DerivedU>
-IGL_INLINE void igl::copyleft::comiso::MIQ_class<DerivedV, DerivedF, DerivedU>::extractUV(Eigen::PlainObjectBase<DerivedU> &UV_out,
-                                                                        Eigen::PlainObjectBase<DerivedF> &FUV_out)
-{
-  UV_out = this->UV_out;
-  FUV_out = this->FUV_out;
-}
-
-template <typename DerivedV, typename DerivedF, typename DerivedU>
-IGL_INLINE int igl::copyleft::comiso::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::copyleft::comiso::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 = Vcut.row(Fcut(f,0));
-  Eigen::Matrix<typename DerivedV::Scalar, 3, 1> p1 = Vcut.row(Fcut(f,1));
-  Eigen::Matrix<typename DerivedV::Scalar, 3, 1> p2 = Vcut.row(Fcut(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::copyleft::comiso::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::copyleft::comiso::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 < Fcut.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);
-
-      stiffnessVector[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::copyleft::comiso::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::copyleft::comiso::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::copyleft::comiso::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::Matrix<int, Eigen::Dynamic, 3> &mismatch,
-  const Eigen::Matrix<int, Eigen::Dynamic, 1> &singular,
-  const Eigen::Matrix<int, Eigen::Dynamic, 3> &seams,
-  Eigen::PlainObjectBase<DerivedU> &UV,
-  Eigen::PlainObjectBase<DerivedF> &FUV,
-  double gradientSize,
-  double stiffness,
-  bool directRound,
-  unsigned int iter,
-  unsigned int localIter,
-  bool doRound,
-  bool singularityRound,
-  const std::vector<int> &roundVertices,
-  const std::vector<std::vector<int>> &hardFeatures)
-{
-  gradientSize = gradientSize/(V.colwise().maxCoeff()-V.colwise().minCoeff()).norm();
-
-  igl::copyleft::comiso::MIQ_class<DerivedV, DerivedF, DerivedU> miq(V,
-    F,
-    PD1_combed,
-    PD2_combed,
-    mismatch,
-    singular,
-    seams,
-    UV,
-    FUV,
-    gradientSize,
-    stiffness,
-    directRound,
-    iter,
-    localIter,
-    doRound,
-    singularityRound,
-    roundVertices,
-    hardFeatures);
-
-  miq.extractUV(UV,FUV);
-}
-
-template <typename DerivedV, typename DerivedF, typename DerivedU>
-IGL_INLINE void igl::copyleft::comiso::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,
-  unsigned int iter,
-  unsigned int localIter,
-  bool doRound,
-  bool singularityRound,
-  const std::vector<int> &roundVertices,
-  const std::vector<std::vector<int>> &hardFeatures)
-{
-
-  DerivedV BIS1, BIS2;
-  igl::compute_frame_field_bisectors(V, F, PD1, PD2, BIS1, BIS2);
-
-  DerivedV BIS1_combed, BIS2_combed;
-  igl::comb_cross_field(V, F, BIS1, BIS2, BIS1_combed, BIS2_combed);
-
-  DerivedF Handle_MMatch;
-  igl::cross_field_mismatch(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, Handle_Seams);
-
-  DerivedV PD1_combed, PD2_combed;
-  igl::comb_frame_field(V, F, PD1, PD2, BIS1_combed, BIS2_combed, PD1_combed, PD2_combed);
-
-  igl::copyleft::comiso::miq(V,
-    F,
-    PD1_combed,
-    PD2_combed,
-    Handle_MMatch,
-    isSingularity,
-    Handle_Seams,
-    UV,
-    FUV,
-    gradientSize,
-    stiffness,
-    directRound,
-    iter,
-    localIter,
-    doRound,
-    singularityRound,
-    roundVertices,
-    hardFeatures);
-}
-
-#ifdef IGL_STATIC_LIBRARY
-// Explicit template instantiation
-template void igl::copyleft::comiso::miq<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const &, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const &, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const &, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const &, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > &, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > &, double, double, bool, unsigned int, unsigned int, bool, bool, const std::vector<int> &, const std::vector<std::vector<int>> &);
-template void igl::copyleft::comiso::miq<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::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::Matrix<int, -1, 3, 0, -1, 3> const &, Eigen::Matrix<int, -1, 1, 0, -1, 1> const &, Eigen::Matrix<int, -1, 3, 0, -1, 3> const &, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > &, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > &, double, double, bool, unsigned int, unsigned int, bool, bool, const std::vector<int> &, const std::vector<std::vector<int>> &);
-template void igl::copyleft::comiso::miq<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::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, -1, 0, -1, -1> > &, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > &, double, double, bool, unsigned int, unsigned int, bool, bool, const std::vector<int> &, const std::vector<std::vector<int>> &);
-#endif

include/igl/copyleft/comiso/miq.h

@@ -1,112 +0,0 @@
-// This file is part of libigl, a simple c++ geometry processing library.
-//
-// Copyright (C) 2014 Daniele Panozzo <[email protected]>, Olga Diamanti <[email protected]>, Kevin Walliman <[email protected]>
-//
-// This Source Code Form is subject to the terms of the Mozilla Public License
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can
-// obtain one at http://mozilla.org/MPL/2.0/.
-#ifndef IGL_COMISO_MIQ_H
-#define IGL_COMISO_MIQ_H
-#include "../../igl_inline.h"
-#include <Eigen/Core>
-#include <vector>
-
-namespace igl
-{
-  namespace copyleft
-  {
-  namespace comiso
-  {
-    /// Global seamless parametrization aligned with a given per-face Jacobian (PD1, PD2).
-    /// The algorithm is based on
-    /// "Mixed-Integer Quadrangulation" by D. Bommes, H. Zimmer, L. Kobbelt
-    /// ACM SIGGRAPH 2009, Article No. 77 (http://dl.acm.org/citation.cfm?id=1531383)
-    /// We thank Nico Pietroni for providing a reference implementation of MIQ
-    /// on which our code is based.
-    ///
-    ///  \bug Due to the way of handling of hardFeatures the algorithm  may fail in difficult cases.
-    ///  \bug Meshes with boundaries are not hendled properly i.e., jagged edges along the boundary are possible
-    /// @param[in] V                 #V by 3 list of mesh vertex 3D positions
-    /// @param[in] F                 #F by 3 list of faces indices in V
-    /// @param[in] PD1               #V by 3 first line of the Jacobian per triangle
-    /// @param[in] PD2               #V by 3 second line of the Jacobian per triangle (optional, if empty it will be a vector in the tangent plane orthogonal to PD1)
-    /// @param[in] gradientSize      global scaling for the gradient (controls the quads resolution)
-    /// @param[in] stiffness         weight for the stiffness iterations (Reserved but not used!)
-    /// @param[in] directRound       greedily round all integer variables at once (greatly improves optimization speed but lowers quality)
-    /// @param[in] iter              stiffness iterations (0 = no stiffness)
-    /// @param[in] localIter         number of local iterations for the integer rounding
-    /// @param[in] doRound           enables the integer rounding (disabling it could be useful for debugging)
-    /// @param[in] singularityRound  set true/false to decide if the singularities' coordinates should be rounded to the nearest integers
-    /// @param[in] roundVertices     id of additional vertices that should be snapped to integer coordinates
-    /// @param[in] hardFeatures      #H by 2 list of pairs of vertices that belongs to edges that should be snapped to integer coordinates
-    /// @param[out] UV                 #UV by 2 list of vertices in 2D
-    /// @param[out] FUV                #FUV by 3 list of face indices in UV
-    ///
-    template <typename DerivedV, typename DerivedF, typename DerivedU>
-    IGL_INLINE void 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 = 30.0,
-      double stiffness = 5.0,
-      bool directRound = false,
-      unsigned int iter = 5,
-      unsigned int localIter = 5,
-      bool doRound = true,
-      bool singularityRound = true,
-      const std::vector<int> &roundVertices = std::vector<int>(),
-      const std::vector<std::vector<int>> &hardFeatures = std::vector<std::vector<int> >());
-
-    /// miq Helper function that allows to directly provided pre-combed bisectors for an already cut mesh
-    ///
-    /// @param[in] V                  #V by 3 list of mesh vertex 3D positions
-    /// @param[in] F                  #F by 3 list of faces indices in V
-    /// @param[in] Additional Input:
-    /// @param[in] PD1_combed         #F by 3 first combed Jacobian
-    /// @param[in] PD2_combed         #F by 3 second combed Jacobian
-    /// @param[in] mismatch             #F by 3 list of per-corner integer PI/2 rotations
-    /// @param[in] singular           #V list of flag that denotes if a vertex is singular or not
-    /// @param[in] seams              #F by 3 list of per-corner flag that denotes seams
-    /// @param[out] UV                 #UV by 2 list of vertices in 2D
-    /// @param[out] FUV                #FUV by 3 list of face indices in UV
-    /// @param[in] gradientSize       global scaling for the gradient (controls the quads resolution)
-    /// @param[in] stiffness          weight for the stiffness iterations (Reserved but not used!)
-    /// @param[in] directRound        greedily round all integer variables at once (greatly improves optimization speed but lowers quality)
-    /// @param[in] iter               stiffness iterations (0 = no stiffness)
-    /// @param[in] localIter          number of local iterations for the integer rounding
-    /// @param[in] doRound            enables the integer rounding (disabling it could be useful for debugging)
-    /// @param[in] singularityRound   set true/false to decide if the singularities' coordinates should be rounded to the nearest integers
-    /// @param[in] roundVertices      id of additional vertices that should be snapped to integer coordinates
-    /// @param[in] hardFeatures       #H by 2 list of pairs of vertices that belongs to edges that should be snapped to integer coordinates
-    ///
-    template <typename DerivedV, typename DerivedF, typename DerivedU>
-    IGL_INLINE void 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::Matrix<int, Eigen::Dynamic, 3> &mismatch,
-      const Eigen::Matrix<int, Eigen::Dynamic, 1> &singular,
-      const Eigen::Matrix<int, Eigen::Dynamic, 3> &seams,
-      Eigen::PlainObjectBase<DerivedU> &UV,
-      Eigen::PlainObjectBase<DerivedF> &FUV,
-      double gradientSize = 30.0,
-      double stiffness = 5.0,
-      bool directRound = false,
-      unsigned int iter = 5,
-      unsigned int localIter = 5,
-      bool doRound = true,
-      bool singularityRound = true,
-      const std::vector<int> &roundVertices = std::vector<int>(),
-      const std::vector<std::vector<int>> &hardFeatures = std::vector<std::vector<int> >());
-  };
-};
-};
-#ifndef IGL_STATIC_LIBRARY
-#include "miq.cpp"
-#endif
-
-#endif

include/igl/copyleft/comiso/nrosy.cpp

@@ -1,784 +0,0 @@
-// This file is part of libigl, a simple c++ geometry processing library.
-//
-// Copyright (C) 2014 Daniele Panozzo <[email protected]>
-//
-// 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 "nrosy.h"
-
-#include "nrosy.h"
-#include "../../triangle_triangle_adjacency.h"
-#include "../../edge_topology.h"
-#include "../../per_face_normals.h"
-
-#include <stdexcept>
-#include "../../PI.h"
-
-#include <Eigen/Geometry>
-#include <Eigen/Sparse>
-#include <queue>
-#include <vector>
-
-#include <gmm/gmm.h>
-#include <CoMISo/Solver/ConstrainedSolver.hh>
-#include <CoMISo/Solver/MISolver.hh>
-#include <CoMISo/Solver/GMM_Tools.hh>
-
-namespace igl
-{
-namespace copyleft
-{
-
-namespace comiso
-{
-class NRosyField
-{
-public:
-  // Init
-  IGL_INLINE NRosyField(const Eigen::MatrixXd& _V, const Eigen::MatrixXi& _F);
-
-  // Generate the N-rosy field
-  // N degree of the rosy field
-  // round separately: round the integer variables one at a time, slower but higher quality
-  IGL_INLINE void solve(int N = 4);
-
-  // Set a hard constraint on fid
-  // fid: face id
-  // v: direction to fix (in 3d)
-  IGL_INLINE void setConstraintHard(int fid, const Eigen::Vector3d& v);
-
-  // Set a soft constraint on fid
-  // fid: face id
-  // w: weight of the soft constraint, clipped between 0 and 1
-  // v: direction to fix (in 3d)
-  IGL_INLINE void setConstraintSoft(int fid, double w, const Eigen::Vector3d& v);
-
-  // 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();
-
-  // Compute singularity indexes
-  IGL_INLINE void findCones(int N);
-
-  // Return the singularities
-  IGL_INLINE Eigen::VectorXd getSingularityIndexPerVertex();
-
-private:
-  // Compute angle differences between reference frames
-  IGL_INLINE void computek();
-
-  // Remove useless matchings
-  IGL_INLINE void reduceSpace();
-
-  // Prepare the system matrix
-  IGL_INLINE void prepareSystemMatrix(int N);
-
-  // Solve with roundings using CoMIso
-  IGL_INLINE void solveRoundings();
-
-  // Convert a vector in 3d to an angle wrt the local reference system
-  IGL_INLINE double convert3DtoLocal(unsigned fid, const Eigen::Vector3d& v);
-
-  // Convert an angle wrt the local reference system to a 3d vector
-  IGL_INLINE Eigen::Vector3d convertLocalto3D(unsigned fid, double a);
-
-  // Compute the per vertex angle defect
-  IGL_INLINE Eigen::VectorXd angleDefect();
-
-  // Temporary variable for the field
-  Eigen::VectorXd angles;
-
-  // Hard constraints
-  Eigen::VectorXd hard;
-  std::vector<bool> isHard;
-
-  // Soft constraints
-  Eigen::VectorXd soft;
-  Eigen::VectorXd wSoft;
-  double softAlpha;
-
-  // Face Topology
-  Eigen::MatrixXi TT, TTi;
-
-  // Edge Topology
-  Eigen::MatrixXi EV, FE, EF;
-  std::vector<bool> isBorderEdge;
-
-  // Per Edge information
-  // Angle between two reference frames
-  Eigen::VectorXd k;
-
-  // Jumps
-  Eigen::VectorXi p;
-  std::vector<bool> pFixed;
-
-  // Mesh
-  Eigen::MatrixXd V;
-  Eigen::MatrixXi F;
-
-  // Normals per face
-  Eigen::MatrixXd N;
-
-  // Singularity index
-  Eigen::VectorXd singularityIndex;
-
-  // Reference frame per triangle
-  std::vector<Eigen::MatrixXd> TPs;
-
-  // System stuff
-  Eigen::SparseMatrix<double> A;
-  Eigen::VectorXd b;
-  Eigen::VectorXi tag_t;
-  Eigen::VectorXi tag_p;
-
-};
-
-} // NAMESPACE COMISO
-} // NAMESPACE COPYLEFT
-} // NAMESPACE IGL
-
-igl::copyleft::comiso::NRosyField::NRosyField(const Eigen::MatrixXd& _V, const Eigen::MatrixXi& _F)
-{
-  V = _V;
-  F = _F;
-
-  assert(V.rows() > 0);
-  assert(F.rows() > 0);
-
-  // Generate topological relations
-  igl::triangle_triangle_adjacency(F,TT,TTi);
-  igl::edge_topology(V,F, EV, FE, EF);
-
-  // Flag border edges
-  isBorderEdge.resize(EV.rows());
-  for(unsigned i=0; i<EV.rows(); ++i)
-    isBorderEdge[i] = (EF(i,0) == -1) || ((EF(i,1) == -1));
-
-  // Generate normals per face
-  igl::per_face_normals(V, F, N);
-
-  // Generate reference frames
-  for(unsigned fid=0; fid<F.rows(); ++fid)
-  {
-    // First edge
-    Eigen::Vector3d e1 = V.row(F(fid,1)) - V.row(F(fid,0));
-    e1.normalize();
-    Eigen::Vector3d e2 = N.row(fid);
-    e2 = e2.cross(e1);
-    e2.normalize();
-
-    Eigen::MatrixXd TP(2,3);
-    TP << e1.transpose(), e2.transpose();
-    TPs.push_back(TP);
-  }
-
-  // Alloc internal variables
-  angles = Eigen::VectorXd::Zero(F.rows());
-  p = Eigen::VectorXi::Zero(EV.rows());
-  pFixed.resize(EV.rows());
-  k = Eigen::VectorXd::Zero(EV.rows());
-  singularityIndex = Eigen::VectorXd::Zero(V.rows());
-
-  // Reset the constraints
-  resetConstraints();
-
-  // Compute k, differences between reference frames
-  computek();
-  softAlpha = 0.5;
-}
-
-void igl::copyleft::comiso::NRosyField::setSoftAlpha(double alpha)
-{
-  assert(alpha >= 0 && alpha < 1);
-  softAlpha = alpha;
-}
-
-
-void igl::copyleft::comiso::NRosyField::prepareSystemMatrix(const int N)
-{
-  double Nd = N;
-
-  // Minimize the MIQ energy
-  // Energy on edge ij is
-  //     (t_i - t_j + kij + pij*(2*pi/N))^2
-  // Partial derivatives:
-  //   t_i: 2     ( t_i - t_j + kij + pij*(2*pi/N)) = 0
-  //   t_j: 2     (-t_i + t_j - kij - pij*(2*pi/N)) = 0
-  //   pij: 4pi/N ( t_i - t_j + kij + pij*(2*pi/N)) = 0
-  //
-  //          t_i      t_j         pij       kij
-  // t_i [     2       -2           4pi/N      2    ]
-  // t_j [    -2        2          -4pi/N     -2    ]
-  // pij [   4pi/N   -4pi/N    2*(2pi/N)^2   4pi/N  ]
-
-  // Count and tag the variables
-  tag_t = Eigen::VectorXi::Constant(F.rows(),-1);
-  std::vector<int> id_t;
-  size_t count = 0;
-  for(unsigned i=0; i<F.rows(); ++i)
-    if (!isHard[i])
-    {
-      tag_t(i) = count++;
-      id_t.push_back(i);
-    }
-
-  size_t count_t = id_t.size();
-
-  tag_p = Eigen::VectorXi::Constant(EF.rows(),-1);
-  std::vector<int> id_p;
-  for(unsigned i=0; i<EF.rows(); ++i)
-  {
-    if (!pFixed[i])
-    {
-      // if it is not fixed then it is a variable
-      tag_p(i) = count++;
-    }
-
-    // if it is not a border edge,
-    if (!isBorderEdge[i])
-    {
-      // and it is not between two fixed faces
-      if (!(isHard[EF(i,0)] && isHard[EF(i,1)]))
-      {
-          // then it participates in the energy!
-          id_p.push_back(i);
-      }
-    }
-  }
-
-  size_t count_p = count - count_t;
-  // System sizes: A (count_t + count_p) x (count_t + count_p)
-  //               b (count_t + count_p)
-
-  b.resize(count_t + count_p);
-  b.setZero();
-
-  std::vector<Eigen::Triplet<double> > T;
-  T.reserve(3 * 4 * count_p);
-
-  for(auto eid : id_p)
-  {
-    int i = EF(eid, 0);
-    int j = EF(eid, 1);
-    bool isFixed_i = isHard[i];
-    bool isFixed_j = isHard[j];
-    bool isFixed_p = pFixed[eid];
-    int row;
-    // (i)-th row: t_i [     2       -2           4pi/N      2    ]
-    if (!isFixed_i)
-    {
-      row = tag_t[i];
-      T.emplace_back(row, tag_t[i], 2);
-      if (isFixed_j)
-        b(row) +=  2 * hard[j];
-      else
-        T.emplace_back(row, tag_t[j], -2);
-      if (isFixed_p)
-        b(row) += -((4. * igl::PI) / Nd) * p[eid];
-      else
-        T.emplace_back(row, tag_p[eid], ((4. * igl::PI) / Nd));
-      b(row) += -2 * k[eid];
-      assert(hard[i] == hard[i]);
-      assert(hard[j] == hard[j]);
-      assert(p[eid] == p[eid]);
-      assert(k[eid] == k[eid]);
-      assert(b(row) == b(row));
-    }
-    // (j)+1 -th row: t_j [    -2        2          -4pi/N     -2    ]
-    if (!isFixed_j)
-    {
-      row = tag_t[j];
-      T.emplace_back(row, tag_t[j], 2);
-      if (isFixed_i)
-        b(row) += 2 * hard[i];
-      else
-        T.emplace_back(row, tag_t[i], -2);
-      if (isFixed_p)
-        b(row) += ((4. * igl::PI) / Nd) * p[eid];
-      else
-        T.emplace_back(row, tag_p[eid], -((4. * igl::PI) / Nd));
-      b(row) += 2 * k[eid];
-      assert(k[eid] == k[eid]);
-      assert(b(row) == b(row));
-    }
-    // (r*3)+2 -th row: pij [   4pi/N   -4pi/N    2*(2pi/N)^2   4pi/N  ]
-    if (!isFixed_p)
-    {
-      row = tag_p[eid];
-      T.emplace_back(row, tag_p[eid], (2. * pow(((2. * igl::PI) / Nd), 2)));
-      if (isFixed_i)
-        b(row) += -(4. * igl::PI) / Nd * hard[i];
-      else
-        T.emplace_back(row, tag_t[i], (4. * igl::PI) / Nd);
-      if (isFixed_j)
-        b(row) += (4. * igl::PI) / Nd * hard[j];
-      else
-        T.emplace_back(row,tag_t[j], -(4. * igl::PI) / Nd);
-      b(row) += - (4 * igl::PI)/Nd * k[eid];
-      assert(k[eid] == k[eid]);
-      assert(b(row) == b(row));
-    }
-  }
-
-  A.resize(count_t + count_p, count_t + count_p);
-  A.setFromTriplets(T.begin(), T.end());
-
-  // Soft constraints
-  bool addSoft = false;
-
-  for(unsigned i=0; i<wSoft.size();++i)
-    if (wSoft[i] != 0)
-      addSoft = true;
-
-  if (addSoft)
-  {
-    Eigen::VectorXd bSoft = Eigen::VectorXd::Zero(count_t + count_p);
-    std::vector<Eigen::Triplet<double> > TSoft;
-    TSoft.reserve(2 * count_p);
-
-    for(unsigned i=0; i<F.rows(); ++i)
-    {
-      int varid = tag_t[i];
-      if (varid != -1) // if it is a variable in the system
-      {
-        TSoft.emplace_back(varid, varid, wSoft[i]);
-        bSoft[varid] += wSoft[i] * soft[i];
-      }
-    }
-    Eigen::SparseMatrix<double> ASoft(count_t + count_p, count_t + count_p);
-    ASoft.setFromTriplets(TSoft.begin(), TSoft.end());
-
-    A = (1.0 - softAlpha) * A + softAlpha * ASoft;
-    b = b * (1.0 - softAlpha) + bSoft * softAlpha;
-  }
-}
-
-void igl::copyleft::comiso::NRosyField::solveRoundings()
-{
-  unsigned n = A.rows();
-
-  gmm::col_matrix< gmm::wsvector< double > > gmm_A(n, n);
-  std::vector<double> gmm_b(n);
-  std::vector<int> ids_to_round;
-  std::vector<double> x(n);
-
-  // Copy A
-  for (int k=0; k<A.outerSize(); ++k)
-    for (Eigen::SparseMatrix<double>::InnerIterator it(A, k); it; ++it)
-    {
-      gmm_A(it.row(),it.col()) += it.value();
-    }
-
-  // Copy b
-  for(unsigned int i = 0; i < n;++i)
-    gmm_b[i] = b[i];
-
-  // Set variables to round
-  ids_to_round.clear();
-  for(unsigned i=0; i<tag_p.size();++i)
-    if(tag_p[i] != -1)
-      ids_to_round.push_back(tag_p[i]);
-
-  // Empty constraints
-  gmm::row_matrix< gmm::wsvector< double > > gmm_C(0, n);
-
-  COMISO::ConstrainedSolver cs;
-  cs.solve(gmm_C, gmm_A, x, gmm_b, ids_to_round, 0.0, false);
-
-
-  // Copy the result back
-  for(unsigned i=0; i<F.rows(); ++i)
-    if (tag_t[i] != -1)
-      angles[i] = x[tag_t[i]];
-    else
-      angles[i] = hard[i];
-
-  for(unsigned i=0; i<EF.rows(); ++i)
-    if(tag_p[i]  != -1)
-      p[i] = (int)std::round(x[tag_p[i]]);
-}
-
-
-void igl::copyleft::comiso::NRosyField::solve(const int N)
-{
-  // Reduce the search space by fixing matchings
-  reduceSpace();
-
-  // Build the system
-  prepareSystemMatrix(N);
-
-  // Solve with integer roundings
-  solveRoundings();
-
-  // Find the cones
-  findCones(N);
-}
-
-void igl::copyleft::comiso::NRosyField::setConstraintHard(const int fid, const Eigen::Vector3d& v)
-{
-  isHard[fid] = true;
-  hard(fid) = convert3DtoLocal(fid, v);
-}
-
-void igl::copyleft::comiso::NRosyField::setConstraintSoft(const int fid, const double w, const Eigen::Vector3d& v)
-{
-  wSoft(fid) = w;
-  soft(fid) = convert3DtoLocal(fid, v);
-}
-
-void igl::copyleft::comiso::NRosyField::resetConstraints()
-{
-  isHard.resize(F.rows());
-  for(unsigned i = 0; i < F.rows(); ++i)
-    isHard[i] = false;
-  hard   = Eigen::VectorXd::Zero(F.rows());
-  wSoft = Eigen::VectorXd::Zero(F.rows());
-  soft = Eigen::VectorXd::Zero(F.rows());
-}
-
-Eigen::MatrixXd igl::copyleft::comiso::NRosyField::getFieldPerFace()
-{
-  Eigen::MatrixXd result(F.rows(),3);
-  for(unsigned int i = 0; i < F.rows(); ++i)
-    result.row(i) = convertLocalto3D(i, angles(i));
-  return result;
-}
-
-void igl::copyleft::comiso::NRosyField::computek()
-{
-  // For every non-border edge
-  for (unsigned eid = 0; eid < EF.rows(); ++eid)
-  {
-    if (!isBorderEdge[eid])
-    {
-      int fid0 = EF(eid,0);
-      int fid1 = EF(eid,1);
-
-      Eigen::Vector3d N0 = N.row(fid0);
-      Eigen::Vector3d N1 = N.row(fid1);
-
-      // find common edge on triangle 0 and 1
-      int fid0_vc = -1;
-      int fid1_vc = -1;
-      for (unsigned i=0;i<3;++i)
-      {
-        if (EV(eid,0) == F(fid0,i))
-          fid0_vc = i;
-        if (EV(eid,1) == F(fid1,i))
-          fid1_vc = i;
-      }
-      assert(fid0_vc != -1);
-      assert(fid1_vc != -1);
-
-      Eigen::Vector3d common_edge = V.row(F(fid0,(fid0_vc+1)%3)) - V.row(F(fid0,fid0_vc));
-      common_edge.normalize();
-
-      // Map the two triangles in a new space where the common edge is the x axis and the N0 the z axis
-      Eigen::MatrixXd P(3,3);
-      Eigen::VectorXd o = V.row(F(fid0,fid0_vc));
-      Eigen::VectorXd tmp = N0.cross(common_edge);
-      P << common_edge, tmp, N0;
-      P.transposeInPlace();
-
-
-      Eigen::MatrixXd V0(3,3);
-      V0.row(0) = V.row(F(fid0,0)).transpose() -o;
-      V0.row(1) = V.row(F(fid0,1)).transpose() -o;
-      V0.row(2) = V.row(F(fid0,2)).transpose() -o;
-
-      V0 = (P*V0.transpose()).transpose();
-
-      assert(V0(0,2) < 1e-10);
-      assert(V0(1,2) < 1e-10);
-      assert(V0(2,2) < 1e-10);
-
-      Eigen::MatrixXd V1(3,3);
-      V1.row(0) = V.row(F(fid1,0)).transpose() -o;
-      V1.row(1) = V.row(F(fid1,1)).transpose() -o;
-      V1.row(2) = V.row(F(fid1,2)).transpose() -o;
-      V1 = (P*V1.transpose()).transpose();
-
-      assert(V1(fid1_vc,2) < 1e-10);
-      assert(V1((fid1_vc+1)%3,2) < 1e-10);
-
-      // compute rotation R such that R * N1 = N0
-      // i.e. map both triangles to the same plane
-      double alpha = -std::atan2(-V1((fid1_vc + 2) % 3, 2), -V1((fid1_vc + 2) % 3, 1));
-
-      Eigen::MatrixXd R(3,3);
-      R << 1,          0,            0,
-           0, std::cos(alpha), -std::sin(alpha),
-           0, std::sin(alpha),  std::cos(alpha);
-      V1 = (R*V1.transpose()).transpose();
-
-      assert(V1(0,2) < 1e-10);
-      assert(V1(1,2) < 1e-10);
-      assert(V1(2,2) < 1e-10);
-
-      // measure the angle between the reference frames
-      // k_ij is the angle between the triangle on the left and the one on the right
-      Eigen::VectorXd ref0 = V0.row(1) - V0.row(0);
-      Eigen::VectorXd ref1 = V1.row(1) - V1.row(0);
-
-      ref0.normalize();
-      ref1.normalize();
-      
-      double ktemp = - std::atan2(ref1(1), ref1(0)) + std::atan2(ref0(1), ref0(0));
-      
-      // make sure kappa is in corret range 
-      auto pos_fmod = [](double x, double y){
-        return (0 == y) ? x : x - y * floor(x/y);
-      };
-      ktemp = pos_fmod(ktemp, 2*igl::PI);
-      if (ktemp > igl::PI) 
-        ktemp -= 2*igl::PI;
-      
-      // just to be sure, rotate ref0 using angle ktemp...
-      Eigen::MatrixXd R2(2,2);
-      R2 << std::cos(-ktemp), -std::sin(-ktemp), std::sin(-ktemp), std::cos(-ktemp);
-
-      tmp = R2*ref0.head<2>();
-
-      assert(tmp(0) - ref1(0) < 1e-10);
-      assert(tmp(1) - ref1(1) < 1e-10);
-
-      k[eid] = ktemp;
-    }
-  }
-
-}
-
-void igl::copyleft::comiso::NRosyField::reduceSpace()
-{
-  // All variables are free in the beginning
-  for(unsigned int i = 0; i < EV.rows(); ++i)
-    pFixed[i] = false;
-
-  std::vector<bool> visited(EV.rows(), false);
-  std::vector<bool> starting(EV.rows(), false);
-
-  std::queue<int> q;
-  for(unsigned int i = 0; i < F.rows(); ++i)
-    if (isHard[i] || wSoft[i] != 0)
-    {
-      q.push(i);
-      starting[i] = true;
-    }
-
-  // Reduce the search space (see MI paper)
-  while (!q.empty())
-  {
-    int c = q.front();
-    q.pop();
-
-    visited[c] = true;
-    for(int i=0; i<3; ++i)
-    {
-      int eid = FE(c,i);
-      int fid = TT(c,i);
-
-      // skip borders
-      if (fid != -1)
-      {
-        assert((EF(eid,0) == c && EF(eid,1) == fid) || (EF(eid,1) == c && EF(eid,0) == fid));
-        // for every neighbouring face
-        if (!visited[fid] && !starting[fid])
-        {
-          pFixed[eid] = true;
-          p[eid] = 0;
-          visited[fid] = true;
-          q.push(fid);
-        }
-      }
-      else
-      {
-        // fix borders
-        pFixed[eid] = true;
-        p[eid] = 0;
-      }
-    }
-  }
-
-  // Force matchings between fixed faces
-  for(unsigned int i = 0; i < F.rows();++i)
-  {
-    if (isHard[i])
-    {
-      for(unsigned int j = 0; j < 3; ++j)
-      {
-        int fid = TT(i,j);
-        if ((fid!=-1) && (isHard[fid]))
-        {
-          // i and fid are adjacent and fixed
-          int eid = FE(i,j);
-          int fid0 = EF(eid,0);
-          int fid1 = EF(eid,1);
-
-          pFixed[eid] = true;
-          p[eid] = (int)std::round(2.0 / igl::PI * (hard(fid1) - hard(fid0) - k(eid)));
-        }
-      }
-    }
-  }
-}
-
-double igl::copyleft::comiso::NRosyField::convert3DtoLocal(unsigned fid, const Eigen::Vector3d& v)
-{
-  // Project onto the tangent plane
-  Eigen::Vector2d vp = TPs[fid] * v;
-
-  // Convert to angle
-  return std::atan2(vp(1), vp(0));
-}
-
-Eigen::Vector3d igl::copyleft::comiso::NRosyField::convertLocalto3D(unsigned fid, double a)
-{
-  Eigen::Vector2d vp(std::cos(a), std::sin(a));
-  return vp.transpose() * TPs[fid];
-}
-
-Eigen::VectorXd igl::copyleft::comiso::NRosyField::angleDefect()
-{
-  Eigen::VectorXd A = Eigen::VectorXd::Constant(V.rows(), 2*igl::PI);
-
-  for (unsigned int i = 0; i < F.rows(); ++i)
-  {
-    for (int j = 0; j < 3; ++j)
-    {
-      Eigen::VectorXd a = V.row(F(i,(j+1)%3)) - V.row(F(i,j));
-      Eigen::VectorXd b = V.row(F(i,(j+2)%3)) - V.row(F(i,j));
-      double t = a.transpose() * b;
-      if(a.norm() > 0. && b.norm() > 0.)
-        t /= (a.norm() * b.norm());
-      else
-        throw std::runtime_error("igl::copyleft::comiso::NRosyField::angleDefect: Division by zero!");
-      A(F(i, j)) -= std::acos(std::max(std::min(t, 1.), -1.));
-    }
-  }
-
-  return A;
-}
-
-void igl::copyleft::comiso::NRosyField::findCones(int N)
-{
-  // Compute I0, see http://www.graphics.rwth-aachen.de/media/papers/bommes_zimmer_2009_siggraph_011.pdf for details
-
-  singularityIndex = Eigen::VectorXd::Zero(V.rows());
-
-  // first the k
-  for (unsigned i = 0; i < EV.rows(); ++i)
-  {
-    if (!isBorderEdge[i])
-    {
-      singularityIndex(EV(i, 0)) += k(i);
-      singularityIndex(EV(i, 1)) -= k(i);
-    }
-  }
-
-  // then the A
-  Eigen::VectorXd A = angleDefect();
-  singularityIndex += A;
-  // normalize
-  singularityIndex /= (2 * igl::PI);
-
-  // round to integer (remove numerical noise)
-  for (unsigned i = 0; i < singularityIndex.size(); ++i)
-    singularityIndex(i) = round(singularityIndex(i));
-
-  for (unsigned i = 0; i < EV.rows(); ++i)
-  {
-    if (!isBorderEdge[i])
-    {
-      singularityIndex(EV(i, 0)) += double(p(i)) / double(N);
-      singularityIndex(EV(i, 1)) -= double(p(i)) / double(N);
-    }
-  }
-
-  // Clear the vertices on the edges
-  for (unsigned i = 0; i < EV.rows(); ++i)
-  {
-    if (isBorderEdge[i])
-    {
-      singularityIndex(EV(i,0)) = 0;
-      singularityIndex(EV(i,1)) = 0;
-    }
-  }
-}
-
-Eigen::VectorXd igl::copyleft::comiso::NRosyField::getSingularityIndexPerVertex()
-{
-  return singularityIndex;
-}
-
-IGL_INLINE void igl::copyleft::comiso::nrosy(
-  const Eigen::MatrixXd& V,
-  const Eigen::MatrixXi& F,
-  const Eigen::VectorXi& b,
-  const Eigen::MatrixXd& bc,
-  const Eigen::VectorXi& b_soft,
-  const Eigen::VectorXd& w_soft,
-  const Eigen::MatrixXd& bc_soft,
-  const int N,
-  const double soft,
-  Eigen::MatrixXd& R,
-  Eigen::VectorXd& S
-  )
-{
-  // Init solver
-  igl::copyleft::comiso::NRosyField solver(V, F);
-
-  // Add hard constraints
-  for (unsigned i = 0; i < b.size(); ++i)
-    solver.setConstraintHard(b(i), bc.row(i));
-
-  // Add soft constraints
-  for (unsigned i = 0; i < b_soft.size(); ++i)
-    solver.setConstraintSoft(b_soft(i), w_soft(i), bc_soft.row(i));
-
-  // Set the soft constraints global weight
-  solver.setSoftAlpha(soft);
-
-  // Interpolate
-  solver.solve(N);
-
-  // Copy the result back
-  R = solver.getFieldPerFace();
-
-  // Extract singularity indices
-  S = solver.getSingularityIndexPerVertex();
-}
-
-
-IGL_INLINE void igl::copyleft::comiso::nrosy(
-                           const Eigen::MatrixXd& V,
-                           const Eigen::MatrixXi& F,
-                           const Eigen::VectorXi& b,
-                           const Eigen::MatrixXd& bc,
-                           const int N,
-                           Eigen::MatrixXd& R,
-                           Eigen::VectorXd& S
-                           )
-{
-  // Init solver
-  igl::copyleft::comiso::NRosyField solver(V, F);
-
-  // Add hard constraints
-  for (unsigned i= 0; i < b.size(); ++i)
-    solver.setConstraintHard(b(i), bc.row(i));
-
-  // Interpolate
-  solver.solve(N);
-
-  // Copy the result back
-  R = solver.getFieldPerFace();
-
-  // Extract singularity indices
-  S = solver.getSingularityIndexPerVertex();
-}

include/igl/copyleft/comiso/nrosy.h

@@ -1,68 +0,0 @@
-// This file is part of libigl, a simple c++ geometry processing library.
-//
-// Copyright (C) 2014 Daniele Panozzo <[email protected]>
-//
-// This Source Code Form is subject to the terms of the Mozilla Public License
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can
-// obtain one at http://mozilla.org/MPL/2.0/.
-#ifndef IGL_COMISO_NROSY_H
-#define IGL_COMISO_NROSY_H
-
-#include <Eigen/Core>
-#include <Eigen/Sparse>
-#include "../../igl_inline.h"
-
-namespace igl
-{
-  namespace copyleft
-  {
-  namespace comiso
-  {
-    /// Generate a N-RoSy field from a sparse set of constraints
-    ///
-    /// @param[in] V       #V by 3 list of mesh vertex coordinates
-    /// @param[in] F       #F by 3 list of mesh faces (must be triangles)
-    /// @param[in] b       #B by 1 list of constrained face indices
-    /// @param[in] bc      #B by 3 list of representative vectors for the
-    ///   constrained faces
-    /// @param[in] b_soft  #S by 1 b for soft constraints
-    /// @param[in] w_soft  #S by 1 weight for the soft constraints (0-1)
-    /// @param[in] bc_soft #S by 3 bc for soft constraints
-    /// @param[in] N       the degree of the N-RoSy vector field
-    /// @param[in] soft    the strength of the soft constraints w.r.t. smoothness
-    ///           (0 -> smoothness only, 1->constraints only)
-    /// @param[out] R       #F by 3 the representative vectors of the interpolated field
-    /// @param[out] S       #V by 1 the singularity index for each vertex (0 = regular)
-    IGL_INLINE void nrosy(
-      const Eigen::MatrixXd& V,
-      const Eigen::MatrixXi& F,
-      const Eigen::VectorXi& b,
-      const Eigen::MatrixXd& bc,
-      const Eigen::VectorXi& b_soft,
-      const Eigen::VectorXd& w_soft,
-      const Eigen::MatrixXd& bc_soft,
-      int N,
-      double soft,
-      Eigen::MatrixXd& R,
-      Eigen::VectorXd& S
-      );
-    /// \overload
-    IGL_INLINE void nrosy(
-     const Eigen::MatrixXd& V,
-     const Eigen::MatrixXi& F,
-     const Eigen::VectorXi& b,
-     const Eigen::MatrixXd& bc,
-     int N,
-     Eigen::MatrixXd& R,
-     Eigen::VectorXd& S
-      );
-
-  }
-}
-}
-
-#ifndef IGL_STATIC_LIBRARY
-#  include "nrosy.cpp"
-#endif
-
-#endif

tests/include/igl/copyleft/comiso/miq.cpp

@@ -1,130 +0,0 @@
-#include <test_common.h>
-#include <igl/avg_edge_length.h>
-#include <igl/barycenter.h>
-#include <igl/comb_cross_field.h>
-#include <igl/comb_frame_field.h>
-#include <igl/compute_frame_field_bisectors.h>
-#include <igl/cross_field_mismatch.h>
-#include <igl/cut_mesh_from_singularities.h>
-#include <igl/find_cross_field_singularities.h>
-#include <igl/local_basis.h>
-#include <igl/readOFF.h>
-#include <igl/rotate_vectors.h>
-#include <igl/copyleft/comiso/miq.h>
-#include <igl/copyleft/comiso/nrosy.h>
-#include <igl/PI.h>
-#include <igl/serialize.h>
-#include <sstream>
-#include <igl/writeDMAT.h>
-
-TEST_CASE("miq: 3_holes", "[igl/copyleft/comiso]")
-{
-using namespace Eigen;
-
-// Input mesh
-Eigen::MatrixXd V;
-Eigen::MatrixXi F;
-
-// Face barycenters
-Eigen::MatrixXd B;
-
-// Cross field
-Eigen::MatrixXd X1,X2;
-
-// Bisector field
-Eigen::MatrixXd BIS1, BIS2;
-
-// Combed bisector
-Eigen::MatrixXd BIS1_combed, BIS2_combed;
-
-// Per-corner, integer mismatches
-Eigen::Matrix<int, Eigen::Dynamic, 3> MMatch;
-
-// Field singularities
-Eigen::Matrix<int, Eigen::Dynamic, 1> isSingularity, singularityIndex;
-
-// Per corner seams
-Eigen::Matrix<int, Eigen::Dynamic, 3> Seams;
-
-// Combed field
-Eigen::MatrixXd X1_combed, X2_combed;
-
-// Global parametrization
-Eigen::MatrixXd UV;
-Eigen::MatrixXi FUV;
-
-// Global parametrization (reference)
-Eigen::MatrixXd UV_ref;
-Eigen::MatrixXi FUV_ref;
-
-// Load a mesh in OFF format
-igl::readOFF(test_common::data_path("3holes.off"), V, F);
-
-double gradient_size = 50;
-double iter = 0;
-double stiffness = 5.0;
-bool direct_round = 0;
-
-// Compute face barycenters
-igl::barycenter(V, F, B);
-
-// Contrain one face
-VectorXi b(1);
-b << 0;
-MatrixXd bc(1, 3);
-bc << 1, 0, 0;
-
-// Create a smooth 4-RoSy field
-VectorXd S;
-igl::copyleft::comiso::nrosy(V, F, b, bc, VectorXi(), VectorXd(), MatrixXd(), 4, 0.5, X1, S);
-
-// Find the orthogonal vector
-MatrixXd B1, B2, B3;
-igl::local_basis(V, F, B1, B2, B3);
-X2 = igl::rotate_vectors(X1, VectorXd::Constant(1, igl::PI / 2), B1, B2);
-
-// Always work on the bisectors, it is more general
-igl::compute_frame_field_bisectors(V, F, X1, X2, BIS1, BIS2);
-
-// Comb the field, implicitly defining the seams
-igl::comb_cross_field(V, F, BIS1, BIS2, BIS1_combed, BIS2_combed);
-
-// Find the integer mismatches
-    igl::cross_field_mismatch(V, F, BIS1_combed, BIS2_combed, true, MMatch);
-
-// Find the singularities
-igl::find_cross_field_singularities(V, F, MMatch, isSingularity, singularityIndex);
-
-// Cut the mesh, duplicating all vertices on the seams
-igl::cut_mesh_from_singularities(V, F, MMatch, Seams);
-
-// Comb the frame-field accordingly
-igl::comb_frame_field(V, F, X1, X2, BIS1_combed, BIS2_combed, X1_combed, X2_combed);
-
-// Global parametrization
-igl::copyleft::comiso::miq(V,
-          F,
-          X1_combed,
-          X2_combed,
-          MMatch,
-          isSingularity,
-          Seams,
-          UV,
-          FUV,
-          gradient_size,
-          stiffness,
-          direct_round,
-          iter,
-          5,
-          true);
-
-  // Refresh the test data
-  // igl::writeDMAT(test_common::data_path("3holes-miq-UV.dmat"),UV);
-  // igl::writeDMAT(test_common::data_path("3holes-miq-FUV.dmat"),FUV);
-
-  igl::readDMAT(test_common::data_path("3holes-miq-UV.dmat"),UV_ref);
-  igl::readDMAT(test_common::data_path("3holes-miq-FUV.dmat"),FUV_ref);
-
-  REQUIRE (1e-6 > (UV-UV_ref).array().abs().maxCoeff());
-  REQUIRE (1e-6 > (FUV-FUV_ref).array().abs().maxCoeff());
-}

tests/include/igl/ray_mesh_intersect.cpp

@@ -30,7 +30,13 @@ TEST_CASE("ray_mesh_intersect: one_triangle", "[igl]")
   IGL_POP_FPE;
 }
 
-TEST_CASE("ray_mesh_intersect: corner-case", "[igl]")
+// https://github.com/libigl/libigl/issues/2363
+#ifdef _WIN32
+#define IGL_MAYFAIL_WIN "[!mayfail]"
+#else
+#define IGL_MAYFAIL_WIN
+#endif
+TEST_CASE("ray_mesh_intersect: corner-case", "[igl]" IGL_MAYFAIL_WIN)
 {
   IGL_PUSH_FPE;
   //       .      //

tutorial/504_NRosyDesign/main.cpp

@@ -1,148 +0,0 @@
-#include <igl/avg_edge_length.h>
-#include <igl/barycenter.h>
-#include <igl/local_basis.h>
-#include <igl/readOFF.h>
-#include <igl/copyleft/comiso/nrosy.h>
-#include <igl/opengl/glfw/Viewer.h>
-#include <igl/PI.h>
-
-
-// Mesh
-Eigen::MatrixXd V;
-Eigen::MatrixXi F;
-
-// Constrained faces id
-Eigen::VectorXi b;
-
-// Cosntrained faces representative vector
-Eigen::MatrixXd bc;
-
-// Degree of the N-RoSy field
-int N = 4;
-
-// Converts a representative vector per face in the full set of vectors that describe
-// an N-RoSy field
-void representative_to_nrosy(
-  const Eigen::MatrixXd& V,
-  const Eigen::MatrixXi& F,
-  const Eigen::MatrixXd& R,
-  const int N,
-  Eigen::MatrixXd& Y)
-{
-  using namespace Eigen;
-  using namespace std;
-  MatrixXd B1, B2, B3;
-
-  igl::local_basis(V,F,B1,B2,B3);
-
-  Y.resize(F.rows()*N,3);
-  for (unsigned i=0;i<F.rows();++i)
-  {
-    double x = R.row(i) * B1.row(i).transpose();
-    double y = R.row(i) * B2.row(i).transpose();
-    double angle = atan2(y,x);
-
-    for (unsigned j=0; j<N;++j)
-    {
-      double anglej = angle + 2*igl::PI*double(j)/double(N);
-      double xj = cos(anglej);
-      double yj = sin(anglej);
-      Y.row(i*N+j) = xj * B1.row(i) + yj * B2.row(i);
-    }
-  }
-}
-
-// Plots the mesh with an N-RoSy field and its singularities on top
-// The constrained faces (b) are colored in red.
-void plot_mesh_nrosy(
-  igl::opengl::glfw::Viewer& viewer,
-  Eigen::MatrixXd& V,
-  Eigen::MatrixXi& F,
-  int N,
-  Eigen::MatrixXd& PD1,
-  Eigen::VectorXd& S,
-  Eigen::VectorXi& b)
-{
-  using namespace Eigen;
-  using namespace std;
-  // Clear the mesh
-  viewer.data().clear();
-  viewer.data().set_mesh(V,F);
-
-  // Expand the representative vectors in the full vector set and plot them as lines
-  double avg = igl::avg_edge_length(V, F);
-  MatrixXd Y;
-  representative_to_nrosy(V, F, PD1, N, Y);
-
-  MatrixXd B;
-  igl::barycenter(V,F,B);
-
-  MatrixXd Be(B.rows()*N,3);
-  for(unsigned i=0; i<B.rows();++i)
-    for(unsigned j=0; j<N; ++j)
-      Be.row(i*N+j) = B.row(i);
-
-  viewer.data().add_edges(Be,Be+Y*(avg/2),RowVector3d(0,0,1));
-
-  // Plot the singularities as colored dots (red for positive, blue for negative)
-  for (unsigned i=0; i<S.size();++i)
-  {
-    if (S(i) < -0.001)
-      viewer.data().add_points(V.row(i),RowVector3d(0,0,1));
-    else if (S(i) > 0.001)
-      viewer.data().add_points(V.row(i),RowVector3d(1,0,0));
-  }
-
-  // Highlight in red the constrained faces
-  MatrixXd C = MatrixXd::Constant(F.rows(),3,1);
-  for (unsigned i=0; i<b.size();++i)
-    C.row(b(i)) << 1, 0, 0;
-  viewer.data().set_colors(C);
-}
-
-  // It allows to change the degree of the field when a number is pressed
-bool key_down(igl::opengl::glfw::Viewer& viewer, unsigned char key, int modifier)
-{
-  using namespace Eigen;
-  using namespace std;
-  if (key >= '1' && key <= '9')
-    N = key - '0';
-
-  MatrixXd R;
-  VectorXd S;
-
-  igl::copyleft::comiso::nrosy(V,F,b,bc,VectorXi(),VectorXd(),MatrixXd(),N,0.5,R,S);
-  plot_mesh_nrosy(viewer,V,F,N,R,S,b);
-
-  return false;
-}
-
-int main(int argc, char *argv[])
-{
-  using namespace std;
-  using namespace Eigen;
-
-  // Load a mesh in OFF format
-  igl::readOFF(TUTORIAL_SHARED_PATH "/bumpy.off", V, F);
-
-  // Threshold faces with high anisotropy
-  b.resize(1);
-  b << 0;
-  bc.resize(1,3);
-  bc << 1,1,1;
-
-  igl::opengl::glfw::Viewer viewer;
-
-  // Interpolate the field and plot
-  key_down(viewer, '4', 0);
-
-  // Plot the mesh
-  viewer.data().set_mesh(V, F);
-  viewer.callback_key_down = &key_down;
-
-  // Disable wireframe
-  viewer.data().show_lines = false;
-
-  // Launch the viewer
-  viewer.launch();
-}

tutorial/505_MIQ/main.cpp

@@ -1,325 +0,0 @@
-#include <igl/avg_edge_length.h>
-#include <igl/barycenter.h>
-#include <igl/comb_cross_field.h>
-#include <igl/comb_frame_field.h>
-#include <igl/compute_frame_field_bisectors.h>
-#include <igl/cross_field_mismatch.h>
-#include <igl/cut_mesh_from_singularities.h>
-#include <igl/find_cross_field_singularities.h>
-#include <igl/local_basis.h>
-#include <igl/readOFF.h>
-#include <igl/rotate_vectors.h>
-#include <igl/copyleft/comiso/miq.h>
-#include <igl/copyleft/comiso/nrosy.h>
-#include <igl/opengl/glfw/Viewer.h>
-#include <igl/PI.h>
-#include <sstream>
-
-#include <igl/serialize.h>
-
-// Input mesh
-Eigen::MatrixXd V;
-Eigen::MatrixXi F;
-
-
-// Face barycenters
-Eigen::MatrixXd B;
-
-// Scale for visualizing the fields
-double global_scale;
-bool extend_arrows = false;
-
-// Cross field
-Eigen::MatrixXd X1,X2;
-
-// Bisector field
-Eigen::MatrixXd BIS1, BIS2;
-
-// Combed bisector
-Eigen::MatrixXd BIS1_combed, BIS2_combed;
-
-// Per-corner, integer mismatches
-Eigen::Matrix<int, Eigen::Dynamic, 3> MMatch;
-
-// Field singularities
-Eigen::Matrix<int, Eigen::Dynamic, 1> isSingularity, singularityIndex;
-
-// Per corner seams
-Eigen::Matrix<int, Eigen::Dynamic, 3> Seams;
-
-// Combed field
-Eigen::MatrixXd X1_combed, X2_combed;
-
-
-// Global parametrization (with seams)
-Eigen::MatrixXd UV_seams;
-Eigen::MatrixXi FUV_seams;
-
-// Global parametrization
-Eigen::MatrixXd UV;
-Eigen::MatrixXi FUV;
-
-
-// Create a texture that hides the integer translation in the parametrization
-void line_texture(Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic> &texture_R,
-                  Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic> &texture_G,
-                  Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic> &texture_B)
-  {
-    unsigned size = 128;
-    unsigned size2 = size/2;
-    unsigned lineWidth = 3;
-    texture_R.setConstant(size, size, 255);
-    for (unsigned i=0; i<size; ++i)
-      for (unsigned j=size2-lineWidth; j<=size2+lineWidth; ++j)
-        texture_R(i,j) = 0;
-    for (unsigned i=size2-lineWidth; i<=size2+lineWidth; ++i)
-      for (unsigned j=0; j<size; ++j)
-        texture_R(i,j) = 0;
-
-    texture_G = texture_R;
-    texture_B = texture_R;
-  }
-
-bool key_down(igl::opengl::glfw::Viewer& viewer, unsigned char key, int modifier)
-{
-  if (key == 'E')
-  {
-    extend_arrows = !extend_arrows;
-  }
-
-  if (key <'1' || key >'8')
-    return false;
-
-  viewer.data().clear();
-  viewer.data().show_lines = false;
-  viewer.data().show_texture = false;
-
-  if (key == '1')
-  {
-    // Cross field
-    viewer.data().set_mesh(V, F);
-    viewer.data().add_edges(extend_arrows ? B - global_scale*X1 : B, B + global_scale*X1 ,Eigen::RowVector3d(1,0,0));
-    viewer.data().add_edges(extend_arrows ? B - global_scale*X2 : B, B + global_scale*X2 ,Eigen::RowVector3d(0,0,1));
-  }
-
-  if (key == '2')
-  {
-    // Bisector field
-    viewer.data().set_mesh(V, F);
-    viewer.data().add_edges(extend_arrows ? B - global_scale*BIS1 : B, B + global_scale*BIS1 ,Eigen::RowVector3d(1,0,0));
-    viewer.data().add_edges(extend_arrows ? B - global_scale*BIS2 : B, B + global_scale*BIS2 ,Eigen::RowVector3d(0,0,1));
-  }
-
-  if (key == '3')
-  {
-    // Bisector field combed
-    viewer.data().set_mesh(V, F);
-    viewer.data().add_edges(extend_arrows ? B - global_scale*BIS1_combed : B, B + global_scale*BIS1_combed ,Eigen::RowVector3d(1,0,0));
-    viewer.data().add_edges(extend_arrows ? B - global_scale*BIS2_combed : B, B + global_scale*BIS2_combed ,Eigen::RowVector3d(0,0,1));
-  }
-
-  if (key == '4')
-  {
-    // Singularities and cuts
-    viewer.data().set_mesh(V, F);
-
-    // Plot cuts
-    int l_count = Seams.sum();
-    Eigen::MatrixXd P1(l_count,3);
-    Eigen::MatrixXd P2(l_count,3);
-
-    for (unsigned i=0; i<Seams.rows(); ++i)
-    {
-      for (unsigned j=0; j<Seams.cols(); ++j)
-      {
-        if (Seams(i,j) != 0)
-        {
-          P1.row(l_count-1) = V.row(F(i,j));
-          P2.row(l_count-1) = V.row(F(i,(j+1)%3));
-          l_count--;
-        }
-      }
-    }
-
-    viewer.data().add_edges(P1, P2, Eigen::RowVector3d(1, 0, 0));
-
-    // Plot the singularities as colored dots (red for negative, blue for positive)
-    for (unsigned i=0; i<singularityIndex.size();++i)
-    {
-      if (singularityIndex(i) < 2 && singularityIndex(i) > 0)
-        viewer.data().add_points(V.row(i),Eigen::RowVector3d(1,0,0));
-      else if (singularityIndex(i) > 2)
-        viewer.data().add_points(V.row(i),Eigen::RowVector3d(0,1,0));
-    }
-
-  }
-
-  if (key == '5')
-  {
-    // Singularities and cuts, original field
-    // Singularities and cuts
-    viewer.data().set_mesh(V, F);
-    viewer.data().add_edges(extend_arrows ? B - global_scale*X1_combed : B, B + global_scale*X1_combed ,Eigen::RowVector3d(1,0,0));
-    viewer.data().add_edges(extend_arrows ? B - global_scale*X2_combed : B, B + global_scale*X2_combed ,Eigen::RowVector3d(0,0,1));
-
-    // Plot cuts
-    int l_count = Seams.sum();
-    Eigen::MatrixXd P1(l_count,3);
-    Eigen::MatrixXd P2(l_count,3);
-
-    for (unsigned i=0; i<Seams.rows(); ++i)
-    {
-      for (unsigned j=0; j<Seams.cols(); ++j)
-      {
-        if (Seams(i,j) != 0)
-        {
-          P1.row(l_count-1) = V.row(F(i,j));
-          P2.row(l_count-1) = V.row(F(i,(j+1)%3));
-          l_count--;
-        }
-      }
-    }
-
-    viewer.data().add_edges(P1, P2, Eigen::RowVector3d(1, 0, 0));
-
-    // Plot the singularities as colored dots (red for negative, blue for positive)
-    for (unsigned i=0; i<singularityIndex.size();++i)
-    {
-      if (singularityIndex(i) < 2 && singularityIndex(i) > 0)
-        viewer.data().add_points(V.row(i),Eigen::RowVector3d(1,0,0));
-      else if (singularityIndex(i) > 2)
-        viewer.data().add_points(V.row(i),Eigen::RowVector3d(0,1,0));
-    }
-  }
-
-  if (key == '6')
-  {
-    // Global parametrization UV
-    viewer.data().set_mesh(UV, FUV);
-    viewer.data().set_uv(UV);
-    viewer.data().show_lines = true;
-  }
-
-  if (key == '7')
-  {
-    // Global parametrization in 3D
-    viewer.data().set_mesh(V, F);
-    viewer.data().set_uv(UV,FUV);
-    viewer.data().show_texture = true;
-  }
-
-  if (key == '8')
-  {
-    // Global parametrization in 3D with seams
-    viewer.data().set_mesh(V, F);
-    viewer.data().set_uv(UV_seams,FUV_seams);
-    viewer.data().show_texture = true;
-  }
-
-  viewer.data().set_colors(Eigen::RowVector3d(1,1,1));
-
-  // Replace the standard texture with an integer shift invariant texture
-  Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic> texture_R, texture_G, texture_B;
-  line_texture(texture_R, texture_G, texture_B);
-  viewer.data().set_texture(texture_R, texture_B, texture_G);
-
-  viewer.core().align_camera_center(viewer.data().V,viewer.data().F);
-
-  return false;
-}
-
-int main(int argc, char *argv[])
-{
-  using namespace Eigen;
-
-  // Load a mesh in OFF format
-  igl::readOFF(TUTORIAL_SHARED_PATH "/3holes.off", V, F);
-
-  double gradient_size = 50;
-  double iter = 0;
-  double stiffness = 5.0;
-  bool direct_round = 0;
-
-  // Compute face barycenters
-  igl::barycenter(V, F, B);
-
-  // Compute scale for visualizing fields
-  global_scale =  .5*igl::avg_edge_length(V, F);
-
-    // Contrain one face
-    VectorXi b(1);
-    b << 0;
-    MatrixXd bc(1, 3);
-    bc << 1, 0, 0;
-
-    // Create a smooth 4-RoSy field
-    VectorXd S;
-    igl::copyleft::comiso::nrosy(V, F, b, bc, VectorXi(), VectorXd(), MatrixXd(), 4, 0.5, X1, S);
-
-    // Find the orthogonal vector
-    MatrixXd B1, B2, B3;
-    igl::local_basis(V, F, B1, B2, B3);
-    X2 = igl::rotate_vectors(X1, VectorXd::Constant(1, igl::PI / 2), B1, B2);
-
-    // Always work on the bisectors, it is more general
-    igl::compute_frame_field_bisectors(V, F, X1, X2, BIS1, BIS2);
-
-    // Comb the field, implicitly defining the seams
-    igl::comb_cross_field(V, F, BIS1, BIS2, BIS1_combed, BIS2_combed);
-
-    // Find the integer mismatches
-    igl::cross_field_mismatch(V, F, BIS1_combed, BIS2_combed, true, MMatch);
-
-    // Find the singularities
-    igl::find_cross_field_singularities(V, F, MMatch, isSingularity, singularityIndex);
-
-    // Cut the mesh, duplicating all vertices on the seams
-    igl::cut_mesh_from_singularities(V, F, MMatch, Seams);
-
-    // Comb the frame-field accordingly
-    igl::comb_frame_field(V, F, X1, X2, BIS1_combed, BIS2_combed, X1_combed, X2_combed);
-
-  // Global parametrization
-  igl::copyleft::comiso::miq(V,
-           F,
-           X1_combed,
-           X2_combed,
-           MMatch,
-           isSingularity,
-           Seams,
-           UV,
-           FUV,
-           gradient_size,
-           stiffness,
-           direct_round,
-           iter,
-           5,
-           true);
-
-// Global parametrization (with seams, only for demonstration)
-igl::copyleft::comiso::miq(V,
-         F,
-         X1_combed,
-         X2_combed,
-         MMatch,
-         isSingularity,
-         Seams,
-         UV_seams,
-         FUV_seams,
-         gradient_size,
-         stiffness,
-         direct_round,
-         iter,
-         5,
-         false);
-
-  // Plot the mesh
-  igl::opengl::glfw::Viewer viewer;
-
-  // Plot the original mesh with a texture parametrization
-  key_down(viewer,'7',0);
-
-  // Launch the viewer
-  viewer.callback_key_down = &key_down;
-  viewer.launch();
-}

tutorial/506_FrameField/main.cpp

@@ -1,255 +0,0 @@
-#include <igl/avg_edge_length.h>
-#include <igl/barycenter.h>
-#include <igl/frame_field_deformer.h>
-#include <igl/frame_to_cross_field.h>
-#include <igl/jet.h>
-#include <igl/local_basis.h>
-#include <igl/readDMAT.h>
-#include <igl/readOBJ.h>
-#include <igl/rotate_vectors.h>
-#include <igl/copyleft/comiso/nrosy.h>
-#include <igl/copyleft/comiso/miq.h>
-#include <igl/copyleft/comiso/frame_field.h>
-#include <igl/opengl/glfw/Viewer.h>
-#include <igl/PI.h>
-
-
-// Input mesh
-Eigen::MatrixXd V;
-Eigen::MatrixXi F;
-
-// Face barycenters
-Eigen::MatrixXd B;
-
-// Scale for visualizing the fields
-double global_scale;
-
-// Input frame field constraints
-Eigen::VectorXi b;
-Eigen::MatrixXd bc1;
-Eigen::MatrixXd bc2;
-
-// Interpolated frame field
-Eigen::MatrixXd FF1, FF2;
-
-// Deformed mesh
-Eigen::MatrixXd V_deformed;
-Eigen::MatrixXd B_deformed;
-
-// Frame field on deformed
-Eigen::MatrixXd FF1_deformed;
-Eigen::MatrixXd FF2_deformed;
-
-// Cross field on deformed
-Eigen::MatrixXd X1_deformed;
-Eigen::MatrixXd X2_deformed;
-
-// Global parametrization
-Eigen::MatrixXd V_uv;
-Eigen::MatrixXi F_uv;
-
-// Create a texture that hides the integer translation in the parametrization
-void line_texture(Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic> &texture_R,
-                  Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic> &texture_G,
-                  Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic> &texture_B)
-{
-  unsigned size = 128;
-  unsigned size2 = size/2;
-  unsigned lineWidth = 3;
-  texture_R.setConstant(size, size, 255);
-  for (unsigned i=0; i<size; ++i)
-    for (unsigned j=size2-lineWidth; j<=size2+lineWidth; ++j)
-      texture_R(i,j) = 0;
-  for (unsigned i=size2-lineWidth; i<=size2+lineWidth; ++i)
-    for (unsigned j=0; j<size; ++j)
-      texture_R(i,j) = 0;
-
-  texture_G = texture_R;
-  texture_B = texture_R;
-}
-
-bool key_down(igl::opengl::glfw::Viewer& viewer, unsigned char key, int modifier)
-{
-  using namespace std;
-  using namespace Eigen;
-
-  if (key <'1' || key >'6')
-    return false;
-
-  viewer.data().clear();
-  viewer.data().show_lines = false;
-  viewer.data().show_texture = false;
-
-  if (key == '1')
-  {
-    // Frame field constraints
-    viewer.data().set_mesh(V, F);
-
-    MatrixXd F1_t = MatrixXd::Zero(FF1.rows(),FF1.cols());
-    MatrixXd F2_t = MatrixXd::Zero(FF2.rows(),FF2.cols());
-    // Highlight in red the constrained faces
-    MatrixXd C = MatrixXd::Constant(F.rows(),3,1);
-    for (unsigned i=0; i<b.size();++i)
-    {
-      C.row(b(i)) << 1, 0, 0;
-      F1_t.row(b(i)) = bc1.row(i);
-      F2_t.row(b(i)) = bc2.row(i);
-    }
-
-    viewer.data().set_colors(C);
-
-    MatrixXd C1,C2;
-    VectorXd K1 = F1_t.rowwise().norm();
-    VectorXd K2 = F2_t.rowwise().norm();
-    igl::jet(K1,true,C1);
-    igl::jet(K2,true,C2);
-
-    viewer.data().add_edges(B - global_scale*F1_t, B + global_scale*F1_t ,C1);
-    viewer.data().add_edges(B - global_scale*F2_t, B + global_scale*F2_t ,C2);
-  }
-
-  if (key == '2')
-  {
-    // Frame field
-    viewer.data().set_mesh(V, F);
-    MatrixXd C1,C2;
-    VectorXd K1 = FF1.rowwise().norm();
-    VectorXd K2 = FF2.rowwise().norm();
-    igl::jet(K1,true,C1);
-    igl::jet(K2,true,C2);
-
-    viewer.data().add_edges(B - global_scale*FF1, B + global_scale*FF1 ,C1);
-    viewer.data().add_edges(B - global_scale*FF2, B + global_scale*FF2 ,C2);
-
-    // Highlight in red the constrained faces
-    MatrixXd C = MatrixXd::Constant(F.rows(),3,1);
-    for (unsigned i=0; i<b.size();++i)
-      C.row(b(i)) << 1, 0, 0;
-    viewer.data().set_colors(C);
-
-  }
-
-  if (key == '3')
-  {
-    // Deformed with frame field
-    viewer.data().set_mesh(V_deformed, F);
-    viewer.data().add_edges(B_deformed - global_scale*FF1_deformed, B_deformed + global_scale*FF1_deformed ,Eigen::RowVector3d(1,0,0));
-    viewer.data().add_edges(B_deformed - global_scale*FF2_deformed, B_deformed + global_scale*FF2_deformed ,Eigen::RowVector3d(0,0,1));
-    viewer.data().set_colors(RowVector3d(1,1,1));
-  }
-
-  if (key == '4')
-  {
-    // Deformed with cross field
-    viewer.data().set_mesh(V_deformed, F);
-    viewer.data().add_edges(B_deformed - global_scale*X1_deformed, B_deformed + global_scale*X1_deformed ,Eigen::RowVector3d(0,0,1));
-    viewer.data().add_edges(B_deformed - global_scale*X2_deformed, B_deformed + global_scale*X2_deformed ,Eigen::RowVector3d(0,0,1));
-    viewer.data().set_colors(RowVector3d(1,1,1));
-  }
-
-  if (key == '5')
-  {
-    // Deformed with quad texture
-    viewer.data().set_mesh(V_deformed, F);
-    viewer.data().set_uv(V_uv,F_uv);
-    viewer.data().set_colors(RowVector3d(1,1,1));
-    viewer.data().show_texture = true;
-  }
-
-  if (key == '6')
-  {
-    // Deformed with quad texture
-    viewer.data().set_mesh(V, F);
-    viewer.data().set_uv(V_uv,F_uv);
-    viewer.data().set_colors(RowVector3d(1,1,1));
-    viewer.data().show_texture = true;
-  }
-
-  // Replace the standard texture with an integer shift invariant texture
-  Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic> texture_R, texture_G, texture_B;
-  line_texture(texture_R, texture_G, texture_B);
-  viewer.data().set_texture(texture_R, texture_B, texture_G);
-  viewer.core().align_camera_center(viewer.data().V,viewer.data().F);
-
-  return false;
-}
-
-int main(int argc, char *argv[])
-{
-  using namespace Eigen;
-
-  // Load a mesh in OBJ format
-  igl::readOBJ(TUTORIAL_SHARED_PATH "/bumpy-cube.obj", V, F);
-
-  // Compute face barycenters
-  igl::barycenter(V, F, B);
-
-  // Compute scale for visualizing fields
-  global_scale =  .2*igl::avg_edge_length(V, F);
-
-  // Load constraints
-  MatrixXd temp;
-  igl::readDMAT(TUTORIAL_SHARED_PATH "/bumpy-cube.dmat",temp);
-
-  b   = temp.block(0,0,temp.rows(),1).cast<int>();
-  bc1 = temp.block(0,1,temp.rows(),3);
-  bc2 = temp.block(0,4,temp.rows(),3);
-
-  // Interpolate the frame field
-  igl::copyleft::comiso::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);
-
-  // Compute face barycenters deformed mesh
-  igl::barycenter(V_deformed, F, B_deformed);
-
-  // Find the closest crossfield to the deformed frame field
-  igl::frame_to_cross_field(V_deformed,F,FF1_deformed,FF2_deformed,X1_deformed);
-
-  // Find a smooth crossfield that interpolates the deformed constraints
-  MatrixXd bc_x(b.size(),3);
-  for (unsigned i=0; i<b.size();++i)
-    bc_x.row(i) = X1_deformed.row(b(i));
-
-  VectorXd S;
-  igl::copyleft::comiso::nrosy(
-             V,
-             F,
-             b,
-             bc_x,
-             VectorXi(),
-             VectorXd(),
-             MatrixXd(),
-             4,
-             0.5,
-             X1_deformed,
-             S);
-
-  // The other representative of the cross field is simply rotated by 90 degrees
-  MatrixXd B1,B2,B3;
-  igl::local_basis(V_deformed,F,B1,B2,B3);
-  X2_deformed =
-    igl::rotate_vectors(X1_deformed, VectorXd::Constant(1,igl::PI/2), B1, B2);
-
-  // Global seamless parametrization
-  igl::copyleft::comiso::miq(V_deformed,
-           F,
-           X1_deformed,
-           X2_deformed,
-           V_uv,
-           F_uv,
-           60.0,
-           5.0,
-           false,
-           2);
-
-  igl::opengl::glfw::Viewer viewer;
-  // Plot the original mesh with a texture parametrization
-  key_down(viewer,'6',0);
-
-  // Launch the viewer
-  viewer.callback_key_down = &key_down;
-  viewer.launch();
-}

tutorial/CMakeLists.txt

@@ -73,10 +73,7 @@ if(LIBIGL_TUTORIALS_CHAPTER5)
     igl_add_tutorial(501_HarmonicParam igl::glfw)
     igl_add_tutorial(502_LSCMParam igl::glfw)
     igl_add_tutorial(503_ARAPParam igl::glfw)
-    igl_add_tutorial(504_NRosyDesign igl::glfw igl_copyleft::comiso)
-    igl_add_tutorial(505_MIQ igl::glfw igl_copyleft::comiso)
-    igl_add_tutorial(506_FrameField igl::glfw igl_copyleft::comiso)
-    igl_add_tutorial(507_Planarization igl::glfw)
+    igl_add_tutorial(504_Planarization igl::glfw)
 endif()
 
 # Chapter 6