diff --git a/python/505_MIQ.py b/python/505_MIQ.py
new file mode 100755
index 000000000..04804fb23
--- /dev/null
+++ b/python/505_MIQ.py
@@ -0,0 +1,279 @@
+import igl
+from math import pi
+
+V = igl.eigen.MatrixXd()
+F = igl.eigen.MatrixXi()
+
+# Face barycenters
+B = igl.eigen.MatrixXd()
+
+# Scale for visualizing the fields
+global_scale = 1
+extend_arrows = False;
+
+# Cross field
+X1 = igl.eigen.MatrixXd()
+X2 = igl.eigen.MatrixXd()
+
+# Bisector field
+BIS1 = igl.eigen.MatrixXd()
+BIS2 = igl.eigen.MatrixXd()
+
+# Combed bisector
+BIS1_combed = igl.eigen.MatrixXd()
+BIS2_combed = igl.eigen.MatrixXd()
+
+# Per-corner, integer mismatches
+MMatch = igl.eigen.MatrixXi()
+
+# Field singularities
+isSingularity = igl.eigen.MatrixXi()
+singularityIndex = igl.eigen.MatrixXi()
+
+# Per corner seams
+Seams = igl.eigen.MatrixXi()
+
+# Combed field
+X1_combed = igl.eigen.MatrixXd()
+X2_combed = igl.eigen.MatrixXd()
+
+# Global parametrization (with seams)
+UV_seams = igl.eigen.MatrixXd()
+FUV_seams = igl.eigen.MatrixXi()
+
+# Global parametrization
+UV = igl.eigen.MatrixXd()
+FUV = igl.eigen.MatrixXi()
+
+# Texture
+texture_R = igl.eigen.MatrixXuc()
+texture_G = igl.eigen.MatrixXuc()
+texture_B = igl.eigen.MatrixXuc()
+
+# Create a texture that hides the integer translation in the parametrization
+def line_texture():
+    size = 128
+    size2 = int(size/2)
+    lineWidth = 3
+    texture_R.setConstant(size, size, 255)
+
+    for i in range(0,size):
+        for j in range(size2-lineWidth,size2+lineWidth+1):
+            texture_R[i,j] = 0
+
+    for i in range(size2-lineWidth,size2+lineWidth+1):
+        for j in range(0,size):
+            texture_R[i,j] = 0
+
+    texture_G = texture_R.copy()
+    texture_B = texture_R.copy()
+    return (texture_R, texture_G, texture_B)
+
+def key_down(viewer, key, modifier):
+    global extend_arrows, texture_R, texture_G, texture_B
+
+    if key == ord('E'):
+        extend_arrows = not extend_arrows
+
+    if key < ord('1') or key > ord('8'):
+        return False;
+
+    viewer.data.clear()
+    viewer.core.show_lines = False
+    viewer.core.show_texture = False
+
+    if key == ord('1'):
+        # Cross field
+        viewer.data.set_mesh(V, F)
+        viewer.data.add_edges(B - global_scale*X1 if extend_arrows else B, B + global_scale*X1 , igl.eigen.MatrixXd([[1,0,0]]))
+        viewer.data.add_edges(B - global_scale*X2 if extend_arrows else B, B + global_scale*X2 , igl.eigen.MatrixXd([[0,0,1]]))
+
+    if key == ord('2'):
+        # Bisector field
+        viewer.data.set_mesh(V, F)
+        viewer.data.add_edges(B - global_scale*BIS1 if extend_arrows else B, B + global_scale*BIS1 , igl.eigen.MatrixXd([[1,0,0]]))
+        viewer.data.add_edges(B - global_scale*BIS2 if extend_arrows else B, B + global_scale*BIS2 , igl.eigen.MatrixXd([[0,0,1]]))
+
+    if key == ord('3'):
+        # Bisector field combed
+        viewer.data.set_mesh(V, F)
+        viewer.data.add_edges(B - global_scale*BIS1_combed if extend_arrows else B, B + global_scale*BIS1_combed , igl.eigen.MatrixXd([[1,0,0]]))
+        viewer.data.add_edges(B - global_scale*BIS2_combed if extend_arrows else B, B + global_scale*BIS2_combed , igl.eigen.MatrixXd([[0,0,1]]))
+
+    if key == ord('4'):
+        # Singularities and cuts
+        viewer.data.set_mesh(V, F)
+
+        # Plot cuts
+        l_count = Seams.sum()
+        P1 = igl.eigen.MatrixXd(l_count,3)
+        P2 = igl.eigen.MatrixXd(l_count,3)
+
+        for i in range(0,Seams.rows()):
+            for j in range(0,Seams.cols()):
+                if Seams[i,j] != 0:
+                    P1.setRow(l_count-1, V.row(F[i,j]))
+                    P2.setRow(l_count-1, V.row(F[i,(j+1)%3]))
+                    l_count = l_count - 1
+
+        viewer.data.add_edges(P1, P2, igl.eigen.MatrixXd([[1, 0, 0]]))
+
+        # Plot the singularities as colored dots (red for negative, blue for positive)
+        for i in range(0,singularityIndex.size()):
+            if singularityIndex[i] < 2 and singularityIndex[i] > 0:
+                viewer.data.add_points(V.row(i),igl.eigen.MatrixXd([[1,0,0]]))
+            elif singularityIndex[i] > 2:
+                viewer.data.add_points(V.row(i),igl.eigen.MatrixXd([[1,0,0]]))
+
+    if key == ord('5'):
+        # Singularities and cuts, original field
+        # Singularities and cuts
+        viewer.data.set_mesh(V, F)
+        viewer.data.add_edges(B - global_scale*X1_combed if extend_arrows else B, B + global_scale*X1_combed ,igl.eigen.MatrixXd([[1,0,0]]))
+        viewer.data.add_edges(B - global_scale*X2_combed if extend_arrows else B, B + global_scale*X2_combed ,igl.eigen.MatrixXd([[0,0,1]]))
+
+        # Plot cuts
+        l_count = Seams.sum()
+
+        P1 = igl.eigen.MatrixXd(l_count,3)
+        P2 = igl.eigen.MatrixXd(l_count,3)
+
+        for i in range(0, Seams.rows()):
+            for j in range(0, Seams.cols()):
+                if Seams[i,j] != 0:
+                    P1.setRow(l_count-1,V.row(F[i,j]))
+                    P2.setRow(l_count-1,V.row(F[i,(j+1)%3]))
+                    l_count = l_count - 1
+
+        viewer.data.add_edges(P1, P2, igl.eigen.MatrixXd([[1, 0, 0]]))
+
+        # Plot the singularities as colored dots (red for negative, blue for positive)
+        for i in range(0,singularityIndex.size()):
+            if singularityIndex[i] < 2 and singularityIndex[i] > 0:
+                viewer.data.add_points(V.row(i),igl.eigen.MatrixXd([[1,0,0]]))
+            elif singularityIndex[i] > 2:
+                viewer.data.add_points(V.row(i),igl.eigen.MatrixXd([[0,1,0]]))
+
+    if key == ord('6'):
+        # Global parametrization UV
+        viewer.data.set_mesh(UV, FUV)
+        viewer.data.set_uv(UV)
+        viewer.core.show_lines = True
+
+    if key == ord('7'):
+        # Global parametrization in 3D
+        viewer.data.set_mesh(V, F)
+        viewer.data.set_uv(UV,FUV)
+        viewer.core.show_texture = True
+
+    if key == ord('8'):
+        # Global parametrization in 3D with seams
+        viewer.data.set_mesh(V, F)
+        viewer.data.set_uv(UV_seams,FUV_seams)
+        viewer.core.show_texture = True
+
+    viewer.data.set_colors(igl.eigen.MatrixXd([[1,1,1]]))
+
+    viewer.data.set_texture(texture_R, texture_B, texture_G)
+
+    viewer.core.align_camera_center(viewer.data.V,viewer.data.F)
+
+    return False
+
+# Load a mesh in OFF format
+igl.readOFF("../tutorial/shared/3holes.off", V, F)
+
+# 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
+b  = igl.eigen.MatrixXi([[0]])
+bc = igl.eigen.MatrixXd([[1,0,0]])
+
+# Create a smooth 4-RoSy field
+S = igl.eigen.MatrixXd()
+
+igl.comiso.nrosy(V,F,b,bc,igl.eigen.MatrixXi(),igl.eigen.MatrixXd(),igl.eigen.MatrixXd(),4,0.5,X1,S)
+
+# Find the the orthogonal vector
+B1 = igl.eigen.MatrixXd()
+B2 = igl.eigen.MatrixXd()
+B3 = igl.eigen.MatrixXd()
+
+igl.local_basis(V,F,B1,B2,B3)
+
+X2 = igl.rotate_vectors(X1, igl.eigen.MatrixXd.Constant(1,1,pi/2), B1, B2)
+
+gradient_size = 50
+iterations = 0
+stiffness = 5.0
+direct_round = False
+
+# 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_missmatch(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.comiso.miq(V,
+           F,
+           X1_combed,
+           X2_combed,
+           MMatch,
+           isSingularity,
+           Seams,
+           UV,
+           FUV,
+           gradient_size,
+           stiffness,
+           direct_round,
+           iterations,
+           5,
+           True,
+           True);
+
+# Global parametrization (with seams, only for demonstration)
+igl.comiso.miq(V,
+         F,
+         X1_combed,
+         X2_combed,
+         MMatch,
+         isSingularity,
+         Seams,
+         UV_seams,
+         FUV_seams,
+         gradient_size,
+         stiffness,
+         direct_round,
+         iterations,
+         5,
+         False);
+
+# Plot the mesh
+viewer = igl.viewer.Viewer()
+
+# Replace the standard texture with an integer shift invariant texture
+(texture_R, texture_G, texture_B) = line_texture()
+
+# Plot the original mesh with a texture parametrization
+key_down(viewer,ord('7'),0)
+
+# Launch the viewer
+viewer.callback_key_down = key_down
+viewer.launch()
diff --git a/python/py_igl/comiso/py_miq.cpp b/python/py_igl/comiso/py_miq.cpp
new file mode 100644
index 000000000..e4981176b
--- /dev/null
+++ b/python/py_igl/comiso/py_miq.cpp
@@ -0,0 +1,59 @@
+m.def("miq", []
+(
+  const Eigen::MatrixXd &V,
+  const Eigen::MatrixXi &F,
+  const Eigen::MatrixXd &PD1,
+  const Eigen::MatrixXd &PD2,
+  Eigen::MatrixXd &UV,
+  Eigen::MatrixXi &FUV,
+  double scale,
+  double stiffness,
+  bool direct_round,
+  int iter,
+  int local_iter,
+  bool DoRound,bool SingularityRound
+  //  std::vector<int> round_vertices,
+  //  std::vector<std::vector<int> > hard_features
+)
+{
+  std::vector<int> round_vertices;
+  std::vector<std::vector<int> > hard_features;
+
+  igl::comiso::miq(V,F,PD1,PD2,UV,FUV,scale,stiffness,direct_round,iter,local_iter,DoRound, SingularityRound, round_vertices, hard_features);
+}, __doc_igl_comiso_miq,
+py::arg("V"), py::arg("F"), py::arg("PD1"), py::arg("PD2"), py::arg("UV"), py::arg("FUV"), py::arg("scale") = 30.0, py::arg("stiffness") = 5.0, py::arg("direct_round") = false, py::arg("iter") = 5, py::arg("local_iter") = 5, py::arg("DoRound") = true, py::arg("SingularityRound") = true
+// , py::arg("round_vertices"), py::arg("hard_features")
+);
+
+m.def("miq", []
+(
+  const Eigen::MatrixXd &V,
+  const Eigen::MatrixXi &F,
+  const Eigen::MatrixXd &PD1_combed,
+  const Eigen::MatrixXd &PD2_combed,
+  const Eigen::MatrixXi &MMatch,
+  const Eigen::MatrixXi &Singular,
+  const Eigen::MatrixXi &Seams,
+  Eigen::MatrixXd &UV,
+  Eigen::MatrixXi &FUV,
+  double GradientSize,
+  double Stiffness,
+  bool DirectRound,
+  int iter,
+  int localIter, bool DoRound, bool SingularityRound
+  // std::vector<int> roundVertices,
+  // std::vector<std::vector<int> > hardFeatures
+)
+{
+  assert_is_VectorX("Singular",Singular);
+
+  std::vector<int> roundVertices;
+  std::vector<std::vector<int> > hardFeatures;
+
+  igl::comiso::miq(V,F,PD1_combed,PD2_combed,MMatch,Singular,Seams,UV,FUV,GradientSize,Stiffness,DirectRound,iter,localIter,DoRound, SingularityRound, roundVertices, hardFeatures);
+}, __doc_igl_comiso_miq,
+py::arg("V"), py::arg("F"), py::arg("PD1_combed"), py::arg("PD2_combed"),
+py::arg("MMatch"), py::arg("Singular"), py::arg("Seams"),
+py::arg("UV"), py::arg("FUV"), py::arg("GradientSize") = 30.0, py::arg("Stiffness") = 5.0, py::arg("DirectRound") = false, py::arg("iter") = 5, py::arg("localIter") = 5, py::arg("DoRound") = true, py::arg("SingularityRound") = true
+// , py::arg("roundVertices"), py::arg("hardFeatures")
+);
diff --git a/python/py_igl/py_comb_cross_field.cpp b/python/py_igl/py_comb_cross_field.cpp
new file mode 100644
index 000000000..8e7a07fb9
--- /dev/null
+++ b/python/py_igl/py_comb_cross_field.cpp
@@ -0,0 +1,13 @@
+m.def("comb_cross_field", []
+(
+  const Eigen::MatrixXd &V,
+  const Eigen::MatrixXi &F,
+  const Eigen::MatrixXd &PD1in,
+  const Eigen::MatrixXd &PD2in,
+  Eigen::MatrixXd &PD1out,
+  Eigen::MatrixXd &PD2out
+)
+{
+  return igl::comb_cross_field(V,F,PD1in,PD2in,PD1out,PD2out);
+}, __doc_igl_comb_cross_field,
+py::arg("V"), py::arg("F"), py::arg("PD1in"), py::arg("PD2in"), py::arg("PD1out"), py::arg("PD2out"));
diff --git a/python/py_igl/py_comb_frame_field.cpp b/python/py_igl/py_comb_frame_field.cpp
new file mode 100644
index 000000000..36e646ac7
--- /dev/null
+++ b/python/py_igl/py_comb_frame_field.cpp
@@ -0,0 +1,15 @@
+m.def("comb_frame_field", []
+(
+  const Eigen::MatrixXd &V,
+  const Eigen::MatrixXi &F,
+  const Eigen::MatrixXd &PD1,
+  const Eigen::MatrixXd &PD2,
+  const Eigen::MatrixXd &BIS1_combed,
+  const Eigen::MatrixXd &BIS2_combed,
+  Eigen::MatrixXd &PD1_combed,
+  Eigen::MatrixXd &PD2_combed
+)
+{
+  return igl::comb_frame_field(V,F,PD1,PD2,BIS1_combed,BIS2_combed,PD1_combed,PD2_combed);
+}, __doc_igl_comb_frame_field,
+py::arg("V"), py::arg("F"), py::arg("PD1"), py::arg("PD2"), py::arg("BIS1_combed"), py::arg("BIS2_combed"), py::arg("PD1_combed"), py::arg("PD2_combed"));
diff --git a/python/py_igl/py_compute_frame_field_bisectors.cpp b/python/py_igl/py_compute_frame_field_bisectors.cpp
new file mode 100644
index 000000000..8cf3c9941
--- /dev/null
+++ b/python/py_igl/py_compute_frame_field_bisectors.cpp
@@ -0,0 +1,29 @@
+m.def("compute_frame_field_bisectors", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  const Eigen::MatrixXd& B1,
+  const Eigen::MatrixXd& B2,
+  const Eigen::MatrixXd& PD1,
+  const Eigen::MatrixXd& PD2,
+  Eigen::MatrixXd& BIS1,
+  Eigen::MatrixXd& BIS2
+)
+{
+  return igl::compute_frame_field_bisectors(V,F,B1,B2,PD1,PD2,BIS1,BIS2);
+}, __doc_igl_compute_frame_field_bisectors,
+py::arg("V"), py::arg("F"), py::arg("B1"), py::arg("B2"), py::arg("PD1"), py::arg("PD2"), py::arg("BIS1"), py::arg("BIS2"));
+
+m.def("compute_frame_field_bisectors", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  const Eigen::MatrixXd& PD1,
+  const Eigen::MatrixXd& PD2,
+  Eigen::MatrixXd& BIS1,
+  Eigen::MatrixXd& BIS2
+)
+{
+  return igl::compute_frame_field_bisectors(V,F,PD1,PD2,BIS1,BIS2);
+}, __doc_igl_compute_frame_field_bisectors,
+py::arg("V"), py::arg("F"), py::arg("PD1"), py::arg("PD2"), py::arg("BIS1"), py::arg("BIS2"));
diff --git a/python/py_igl/py_cross_field_missmatch.cpp b/python/py_igl/py_cross_field_missmatch.cpp
new file mode 100644
index 000000000..25b85f4ab
--- /dev/null
+++ b/python/py_igl/py_cross_field_missmatch.cpp
@@ -0,0 +1,13 @@
+m.def("cross_field_missmatch", []
+(
+  const Eigen::MatrixXd &V,
+  const Eigen::MatrixXi &F,
+  const Eigen::MatrixXd &PD1,
+  const Eigen::MatrixXd &PD2,
+  const bool isCombed,
+  Eigen::MatrixXi &missmatch
+)
+{
+  return igl::cross_field_missmatch(V,F,PD1,PD2,isCombed,missmatch);
+}, __doc_igl_cross_field_missmatch,
+py::arg("V"), py::arg("F"), py::arg("PD1"), py::arg("PD2"), py::arg("isCombed"), py::arg("missmatch"));
diff --git a/python/py_igl/py_cut_mesh_from_singularities.cpp b/python/py_igl/py_cut_mesh_from_singularities.cpp
new file mode 100644
index 000000000..7ab5b71c3
--- /dev/null
+++ b/python/py_igl/py_cut_mesh_from_singularities.cpp
@@ -0,0 +1,11 @@
+m.def("cut_mesh_from_singularities", []
+(
+  const Eigen::MatrixXd &V,
+  const Eigen::MatrixXi &F,
+  const Eigen::MatrixXi &MMatch,
+  Eigen::MatrixXi &seams
+)
+{
+  return igl::cut_mesh_from_singularities(V,F,MMatch,seams);
+}, __doc_igl_cut_mesh_from_singularities,
+py::arg("V"), py::arg("F"), py::arg("MMatch"), py::arg("seams"));
diff --git a/python/py_igl/py_find_cross_field_singularities.cpp b/python/py_igl/py_find_cross_field_singularities.cpp
new file mode 100644
index 000000000..e7dc5134a
--- /dev/null
+++ b/python/py_igl/py_find_cross_field_singularities.cpp
@@ -0,0 +1,27 @@
+m.def("find_cross_field_singularities", []
+(
+  const Eigen::MatrixXd &V,
+  const Eigen::MatrixXi &F,
+  const Eigen::MatrixXi &Handle_MMatch,
+  Eigen::MatrixXi &isSingularity,
+  Eigen::MatrixXi &singularityIndex
+)
+{
+  return igl::find_cross_field_singularities(V,F,Handle_MMatch,isSingularity,singularityIndex);
+}, __doc_igl_find_cross_field_singularities,
+py::arg("V"), py::arg("F"), py::arg("Handle_MMatch"), py::arg("isSingularity"), py::arg("singularityIndex"));
+
+m.def("find_cross_field_singularities", []
+(
+  const Eigen::MatrixXd &V,
+  const Eigen::MatrixXi &F,
+  const Eigen::MatrixXd &PD1,
+  const Eigen::MatrixXd &PD2,
+  Eigen::MatrixXi &isSingularity,
+  Eigen::MatrixXi &singularityIndex,
+  bool isCombed
+)
+{
+  return igl::find_cross_field_singularities(V,F,PD1,PD2,isSingularity,singularityIndex,isCombed);
+}, __doc_igl_find_cross_field_singularities,
+py::arg("V"), py::arg("F"), py::arg("PD1"), py::arg("PD2"), py::arg("isSingularity"), py::arg("singularityIndex"),  py::arg("isCombed") = false);
diff --git a/python/py_igl/py_rotate_vectors.cpp b/python/py_igl/py_rotate_vectors.cpp
new file mode 100644
index 000000000..81a6eac7f
--- /dev/null
+++ b/python/py_igl/py_rotate_vectors.cpp
@@ -0,0 +1,12 @@
+m.def("rotate_vectors", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXd& A,
+  const Eigen::MatrixXd& B1,
+  const Eigen::MatrixXd& B2
+)
+{
+  assert_is_VectorX("A",A);
+  return igl::rotate_vectors(V,A,B1,B2);
+}, __doc_igl_rotate_vectors,
+py::arg("V"), py::arg("A"), py::arg("B1"), py::arg("B2"));
diff --git a/python/py_vector.cpp b/python/py_vector.cpp
index c92b6eb92..08a3fbd42 100644
--- a/python/py_vector.cpp
+++ b/python/py_vector.cpp
@@ -284,6 +284,9 @@ py::class_<Type> bind_eigen_2(py::module &m, const char *name,
         {
           return Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic>(Eigen::Map<const Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic>>(m.data(),r,c));
         })
+
+        .def("copy", [](const Type &m) { return Type(m); })
+
         ;
     return matrix;
 }