diff --git a/python/504_NRosyDesign.py b/python/504_NRosyDesign.py
new file mode 100755
index 000000000..08f6a4862
--- /dev/null
+++ b/python/504_NRosyDesign.py
@@ -0,0 +1,108 @@
+import igl
+from math import atan2,pi,cos,sin
+
+# Mesh
+V = igl.eigen.MatrixXd()
+F = igl.eigen.MatrixXi()
+
+# Constrained faces id
+b = igl.eigen.MatrixXi()
+
+# Constrained faces representative vector
+bc = igl.eigen.MatrixXd()
+
+# Degree of the N-RoSy field
+N = 4;
+
+# Converts a representative vector per face in the full set of vectors that describe
+# an N-RoSy field
+def representative_to_nrosy(V, F, R, N, Y):
+    B1 = igl.eigen.MatrixXd()
+    B2 = igl.eigen.MatrixXd()
+    B3 = igl.eigen.MatrixXd()
+
+    igl.local_basis(V,F,B1,B2,B3)
+
+    Y.resize(F.rows()*N, 3)
+
+    for i in range (0,F.rows()):
+        x = R.row(i) * B1.row(i).transpose()
+        y = R.row(i) * B2.row(i).transpose()
+        angle = atan2(y[0],x[0])
+
+        for j in range(0,N):
+            anglej = angle + 2*pi*j/float(N)
+            xj = cos(anglej)
+            yj = sin(anglej)
+            Y.setRow(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.
+def plot_mesh_nrosy(viewer, V, F, N, PD1, S, b):
+    # 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
+    avg = igl.avg_edge_length(V, F)
+    Y = igl.eigen.MatrixXd()
+    representative_to_nrosy(V, F, PD1, N, Y)
+
+    B = igl.eigen.MatrixXd()
+    igl.barycenter(V,F,B)
+
+    Be = igl.eigen.MatrixXd(B.rows()*N,3)
+    for i in range(0,B.rows()):
+        for j in range(0,N):
+            Be.setRow(i*N+j,B.row(i))
+
+    viewer.data.add_edges(Be,Be+Y*(avg/2),igl.eigen.MatrixXd([[0,0,1]]))
+
+    # Plot the singularities as colored dots (red for negative, blue for positive)
+    for i in range(0,S.size()):
+        if S[i] < -0.001:
+            viewer.data.add_points(V.row(i),igl.eigen.MatrixXd([[1,0,0]]))
+        elif S[i] > 0.001:
+            viewer.data.add_points(V.row(i),igl.eigen.MatrixXd([[0,1,0]]));
+
+    # Highlight in red the constrained faces
+    C = igl.eigen.MatrixXd.Constant(F.rows(),3,1)
+    for i in range(0,b.size()):
+        C.setRow(b[i], igl.eigen.MatrixXd([[1, 0, 0]]))
+    viewer.data.set_colors(C)
+
+# It allows to change the degree of the field when a number is pressed
+def key_down(viewer, key, modifier):
+    global N
+    if key >= ord('1') and key <= ord('9'):
+        N = key - '0'
+
+    R = igl.eigen.MatrixXd()
+    S = igl.eigen.MatrixXd()
+
+    igl.comiso.nrosy(V,F,b,bc,igl.eigen.MatrixXi(),igl.eigen.MatrixXd(),igl.eigen.MatrixXd(),N,0.5,R,S)
+    plot_mesh_nrosy(viewer,V,F,N,R,S,b)
+
+    return False
+
+# Load a mesh in OFF format
+igl.readOFF("../tutorial/shared/bumpy.off", V, F);
+
+# Threshold faces with high anisotropy
+b  = igl.eigen.MatrixXi([[0]])
+bc = igl.eigen.MatrixXd([[1,1,1]])
+
+viewer = igl.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.core.show_lines = False
+
+# Launch the viewer
+viewer.launch()
diff --git a/python/py_igl.cpp b/python/py_igl.cpp
index 619e30d9c..60a6d9241 100644
--- a/python/py_igl.cpp
+++ b/python/py_igl.cpp
@@ -39,6 +39,7 @@
 #include <igl/boundary_loop.h>
 #include <igl/map_vertices_to_circle.h>
 #include <igl/lscm.h>
+#include <igl/local_basis.h>
 
 void python_export_igl(py::module &m)
 {
@@ -79,4 +80,5 @@ void python_export_igl(py::module &m)
 #include "py_igl/py_boundary_loop.cpp"
 #include "py_igl/py_map_vertices_to_circle.cpp"
 #include "py_igl/py_lscm.cpp"
+#include "py_igl/py_local_basis.cpp"
 }
diff --git a/python/py_igl/comiso/py_nrosy.cpp b/python/py_igl/comiso/py_nrosy.cpp
new file mode 100644
index 000000000..60a0db8aa
--- /dev/null
+++ b/python/py_igl/comiso/py_nrosy.cpp
@@ -0,0 +1,41 @@
+m.def("nrosy", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  const Eigen::MatrixXi& b,
+  const Eigen::MatrixXd& bc,
+  const Eigen::MatrixXi& b_soft,
+  const Eigen::MatrixXd& w_soft,
+  const Eigen::MatrixXd& bc_soft,
+  const int N,
+  const double soft,
+  Eigen::MatrixXd& R,
+  Eigen::MatrixXd& S
+)
+{
+  assert_is_VectorX("b",b);
+  assert_is_VectorX("b_soft",b_soft);
+  assert_is_VectorX("w_soft",w_soft);
+  Eigen::VectorXd St;
+  igl::comiso::nrosy(V,F,b,bc,b_soft,w_soft,bc_soft,N,soft,R,St);
+  S = St;
+}, __doc_igl_comiso_nrosy,
+py::arg("V"), py::arg("F"), py::arg("b"), py::arg("bc"), py::arg("b_soft"), py::arg("w_soft"), py::arg("bc_soft"), py::arg("N"), py::arg("soft"), py::arg("R"), py::arg("S"));
+
+m.def("nrosy", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  const Eigen::MatrixXi& b,
+  const Eigen::MatrixXd& bc,
+  const int N,
+  Eigen::MatrixXd& R,
+  Eigen::MatrixXd& S
+)
+{
+  assert_is_VectorX("b",b);
+  Eigen::VectorXd St;
+  igl::comiso::nrosy(V,F,b,bc,N,R,St);
+  S = St;
+}, __doc_igl_comiso_nrosy,
+py::arg("V"), py::arg("F"), py::arg("b"), py::arg("bc"), py::arg("N"), py::arg("R"), py::arg("S"));
diff --git a/python/py_igl/py_local_basis.cpp b/python/py_igl/py_local_basis.cpp
new file mode 100644
index 000000000..41aa481e2
--- /dev/null
+++ b/python/py_igl/py_local_basis.cpp
@@ -0,0 +1,12 @@
+m.def("local_basis", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  Eigen::MatrixXd& B1,
+  Eigen::MatrixXd& B2,
+  Eigen::MatrixXd& B3
+)
+{
+  return igl::local_basis(V,F,B1,B2,B3);
+}, __doc_igl_local_basis,
+py::arg("V"), py::arg("F"), py::arg("B1"), py::arg("B2"), py::arg("B3"));