diff --git a/README.md b/README.md
index 5d38be1..9b3588b 100644
--- a/README.md
+++ b/README.md
@@ -25,11 +25,13 @@ Python:
 * networkx
 * matplotlib
 * openmesh
+* libigl Python bindings
 * meshplot: only for jupyter notebook 
 
 Install with 
 ```shell
 pip install numpy scipy networkx matplotlib openmesh
+conda install -c conda-forge igl
 conda install -c conda-forge meshplot
 ```
 
diff --git a/meshutility/mesh_cut.py b/meshutility/mesh_cut.py
index 6d55a48..9a7d92b 100644
--- a/meshutility/mesh_cut.py
+++ b/meshutility/mesh_cut.py
@@ -12,13 +12,14 @@
 import numpy as np
 import openmesh as om
 import networkx as nx
+import igl
 
 from .mesh2graph import ff_graph
 
 __all__ = ['cut_along_curve']
 
 
-def cut_along_curve(V, F, curve_index):
+def cut_along_curve(V, F, curve_index, method='igl'):
     """
     Cut mesh along curve.
 
@@ -26,13 +27,43 @@ def cut_along_curve(V, F, curve_index):
     ----------
     V : numpy.array
     F : numpy.array
-    curve_index : list
+    curve_index : list (int)
     Returns
     -------
     mesh : openmesh.TriMesh
     index : list (int)
     """
+    if method == 'igl':
+        return cut_along_curve_igl(V, F, curve_index)
+    else:
+        return cut_along_curve_mu(V, F, curve_index)
+
 
+def cut_along_curve_igl(V, F, curve_index):
+    mesh = om.TriMesh(V, F)
+    cuts = np.zeros_like(F, 'i')
+    vhs = [mesh.vertex_handle(i) for i in curve_index]
+    for i in range(len(curve_index)-1):
+        he = mesh.find_halfedge(vhs[i], vhs[i+1])
+        for _ in range(2):  # two halfedges
+            if not he.is_valid():
+                break
+            if not mesh.is_boundary(he):
+                fh = mesh.face_handle(he)
+                f = F[fh.idx()]
+                for j in range(3):
+                    if {f[j], f[(j+1)%3]} == {curve_index[i], curve_index[i+1]}:  # compare two sets
+                        cuts[fh.idx(),j] = 1
+                        break
+            he = mesh.opposite_halfedge_handle(he)
+
+    V1, F1 = igl.cut_mesh(V, F, cuts)
+    mesh = om.TriMesh(V1, F1)
+    curve_new_index = [F1[np.where(F==i)].max() for i in curve_index]
+    return mesh, curve_new_index
+
+
+def cut_along_curve_mu(V, F, curve_index):
     mesh = om.TriMesh(V, F)
     assert mesh.n_vertices() == V.shape[0], "Invalid input mesh"
     assert(len(curve_index) > 1)
diff --git a/tests/test_split_cut.py b/tests/test_split_cut.py
index aa0eab5..7aedf76 100644
--- a/tests/test_split_cut.py
+++ b/tests/test_split_cut.py
@@ -6,7 +6,7 @@
 
 def test_cut():
     mesh = om.read_trimesh('../data/sphere_cvt.obj')
-    field = mesh.points()[:, 2]
+    field = mesh.points()[:, 1] - mesh.points()[:, 2]
 
     isocurve = mu.IsoCurve(mesh.points(), mesh.face_vertex_indices(), field)
     pts, on_edges, ratios, isocurve_indices = isocurve.extract(0.)
@@ -109,5 +109,5 @@ def intersect(A,B,C,D):
     
 
 if __name__ == '__main__':
-    # test_cut()
-    test_split_rect_1()
+    test_cut()
+    # test_split_rect_1()