diff --git a/src/gpytoolbox/__init__.py b/src/gpytoolbox/__init__.py
index d27209f0..0a651550 100644
--- a/src/gpytoolbox/__init__.py
+++ b/src/gpytoolbox/__init__.py
@@ -109,7 +109,9 @@
 from .linear_blend_skinning import linear_blend_skinning
 from .barycenters import barycenters
 from .sdf_flow import sdf_flow
+from .cut_edges import cut_edges
 from .biharmonic_energy import biharmonic_energy
 from .biharmonic_energy_intrinsic import biharmonic_energy_intrinsic
 from .adjacency_matrix import adjacency_matrix
+from .non_manifold_edges import non_manifold_edges
 from .connected_components import connected_components
diff --git a/src/gpytoolbox/boundary_loops.py b/src/gpytoolbox/boundary_loops.py
index af19cd67..20750d8a 100644
--- a/src/gpytoolbox/boundary_loops.py
+++ b/src/gpytoolbox/boundary_loops.py
@@ -1,13 +1,10 @@
 import numpy as np
 from gpytoolbox.boundary_edges import boundary_edges
+from gpytoolbox.non_manifold_edges import non_manifold_edges
 
 
 def boundary_loops(f, allow_wrong_orientations=True):
-    """Computes oriented boundary loop for each boundary component of a triangle
-    mesh.
-    This function only works on manifold triangle meshes.
-
-    TODO: assert manifoldness when running this function
+    """Computes a list containing the oriented boundary loop for each boundary component of a triangle mesh in the style of a sorted polyline. This function only works on connected (i.e., single component) manifold triangle meshes.
 
     Parameters
     ----------
@@ -33,6 +30,9 @@ def boundary_loops(f, allow_wrong_orientations=True):
     assert f.shape[0] > 0
     assert f.shape[1] == 3
 
+    # check mesh is manifold
+    assert len(non_manifold_edges(f)) == 0, "Mesh is not manifold"
+
     bE = boundary_edges(f)
 
     #Loop through each boundary, edge, marking them as seen, until all have
diff --git a/src/gpytoolbox/cut_edges.py b/src/gpytoolbox/cut_edges.py
new file mode 100644
index 00000000..dae71255
--- /dev/null
+++ b/src/gpytoolbox/cut_edges.py
@@ -0,0 +1,105 @@
+import numpy as np
+import scipy as sp
+
+from .halfedges import halfedges
+from .array_correspondence import array_correspondence
+from .remove_unreferenced import remove_unreferenced
+
+def cut_edges(F,E):
+    """Cut a triangle mesh along a specified set of edges.
+
+    Given a triangle mesh and a set of edges, this returns a new mesh that has been "cut" along those edges; meaning, such that the two faces incident on that edge are no longer combinatorially connected. This is done by duplicating vertices along the cut edges (see note), and creating a new geometrically identical edge between the duplicated vertices.
+
+    Parameters
+    ----------
+    F : (m,3) numpy int array
+        face index list of a triangle mesh
+    E : (k,2) numpy int array
+        index list of edges to cut, indexing the same array as F.
+        If E contains edges that are not present in F, they will be ignored.
+
+    Returns
+    -------
+    G : (m,3) numpy int array
+        face index list of cut triangle mesh
+    I : (l,) numpy int array
+        list of indices into V of vertices in new mesh such that V[I,:] are the
+        vertices in the new mesh.
+        This takes care of correctly duplicating vertices.
+
+    Notes
+    -----
+    Since only vertices that no longer share an edge in common are duplicated, you cannot cut a single interior edge. This function mirrors gptoolbox's cut_edges (https://github.com/alecjacobson/gptoolbox/blob/master/mesh/cut_edges.m)
+
+    Examples
+    --------
+    ```python
+    _,F = gpy.read_mesh("mesh.obj")
+    E = np.array([[0,1], [1,2]])
+    G,I = gpy.cut_edges(F,G)
+    W = V[I,:]
+    # My new mesh is now W,G
+    ```
+    """
+
+    assert F.shape[0]>0, "F must be nonempty."
+    assert F.shape[1]==3, "Only works for triangle meshes."
+    n = np.max(F)+1
+    if E.size==0:
+        return np.arange(F.size[0]), np.arange(n)
+    assert E.shape[1]==2, "E is a (k,2) matrix."
+
+    # This code is a translation of https://github.com/alecjacobson/gptoolbox/blob/master/mesh/cut_edges.m by Alec Jacobson
+    he = halfedges(F)
+    flat_he = np.concatenate([he[:,0,:],he[:,1,:],he[:,2,:]], axis=0)
+    sorted_he = np.sort(flat_he, axis=1)
+    unique_he, unique_inverse = np.unique(sorted_he, axis=0,
+        return_inverse=True)
+    unique_he_to_F = sp.sparse.csr_matrix(
+        (np.ones(unique_inverse.shape[0]),
+            (unique_inverse,
+                np.tile(np.arange(F.shape[0]),3))),
+        shape=(unique_he.shape[0], F.shape[0])
+        )
+
+    FF = np.arange(3*F.shape[0]).reshape((-1,3), order='F')
+    sorted_unique_he = np.sort(unique_he, axis=1)
+    sorted_E = np.sort(E, axis=1)
+    isin_unique_he_but_not_E = np.nonzero(
+        array_correspondence(sorted_unique_he, sorted_E, axis=0) < 0)[0]
+    noncut = sp.sparse.csr_matrix(
+        (np.ones(isin_unique_he_but_not_E.shape[0]),
+            (isin_unique_he_but_not_E, isin_unique_he_but_not_E)),
+        shape=(unique_he.shape[0], unique_he.shape[0])
+        )
+    unique_he_to_EE = sp.sparse.csr_matrix(
+        (np.ones(unique_inverse.shape[0]),
+            (unique_inverse, np.arange(3*F.shape[0]))),
+        shape=(unique_he.shape[0], 3*F.shape[0])
+        )
+    I = np.arange(3*F.shape[0]).reshape(F.shape[0], 3, order='F')
+    VV_to_EE = sp.sparse.csr_matrix(
+        (np.ones(6*F.shape[0]),
+            (np.concatenate((FF[:,0],FF[:,1],FF[:,2],FF[:,0],FF[:,1],FF[:,2])),
+                np.concatenate((I[:,1],I[:,2],I[:,0],I[:,2],I[:,0],I[:,1])))),
+        shape=(3*F.shape[0], 3*F.shape[0])
+        )
+    VV_to_V = sp.sparse.csr_matrix(
+        (np.ones(I.size), (I.ravel(), F.ravel())),
+        shape=(3*F.shape[0], n)
+        )
+    A = (VV_to_EE * (unique_he_to_EE.T * noncut * unique_he_to_EE)
+        * VV_to_EE.T).multiply(VV_to_V * VV_to_V.T)
+
+    I = F.flatten(order='F')
+    VV = np.zeros(np.max(FF)+1) #dummy vertex data for remove unreferenced
+    _,labels = sp.sparse.csgraph.connected_components(A, return_labels=True)
+    VV[labels] = labels
+    I[labels] = I
+    FF = labels[FF]
+    W,_,IM,_ = remove_unreferenced(VV[:,None], FF, return_maps=True)
+    I = I[IM.ravel()<W.shape[0]]
+    G = IM.ravel(order='F')[FF]
+
+    return G,I
+
diff --git a/src/gpytoolbox/marching_squares.py b/src/gpytoolbox/marching_squares.py
index e99aacb2..74da80a7 100644
--- a/src/gpytoolbox/marching_squares.py
+++ b/src/gpytoolbox/marching_squares.py
@@ -3,7 +3,7 @@
 
 def marching_squares(S,GV,nx,ny):
     """
-    Marching squares algorithm for extracting isocontours from a scalar field.
+    Marching squares algorithm for extracting isocontours from a scalar field. Output is given as a list of (unordered) vertices and edge indices into the vertex list.
 
     Parameters
     ----------
@@ -45,7 +45,11 @@ def marching_squares(S,GV,nx,ny):
     # Plot
     plt.figure()
     plt.imshow(S.reshape((nx,ny),order='F'),extent=[-1,1,-1,1])
-    plt.plot(V[:,0],V[:,1],'k')
+    for i in range(E.shape[0]):
+            plt.plot([V[E[i,0],0],V[E[i,1],0]],
+                     [V[E[i,0],1],V[E[i,1],1]],
+                     'k-')
+    plt.show()
     plt.axis('equal')
     plt.show()
     ```
diff --git a/src/gpytoolbox/non_manifold_edges.py b/src/gpytoolbox/non_manifold_edges.py
new file mode 100644
index 00000000..7da15235
--- /dev/null
+++ b/src/gpytoolbox/non_manifold_edges.py
@@ -0,0 +1,43 @@
+import numpy as np
+from gpytoolbox.halfedges import halfedges
+
+def non_manifold_edges(F):
+    """Given a triangle mesh with face indices F, returns (unoriented) edges that are non-manifold; i.e., edges that are incident to more than two faces.
+
+    Parameters
+    ----------
+    F : (m,3) numpy int array
+        face index list of a triangle mesh
+
+    Returns
+    -------
+    ne : (n,2) numpy int array
+        list of unique non-manifold edges. Columns are sorted in ascending index order, and rows are sorted in lexicographic order.
+
+    Notes
+    -----
+    It would be nice to also have a non_manifold_vertices function that wraps 2D and 3D functionality.
+
+    Examples
+    --------
+    ```python
+    from gpy import non_manifold_edges
+    # bad mesh with one non-manifold edge in [0,2]
+    f = np.array([[0,1,2],[0,2,3],[2,0,4]],dtype=int)
+    ne = gpy.non_manifold_edges(f)
+    # ne is now np.array([[0,2]])
+    ```
+    
+    """
+
+    assert F.shape[1] == 3
+
+    he = halfedges(F).reshape(-1,2)
+    he = np.sort(he, axis=1)
+    # print(he)
+    he_u = np.unique(he, axis=0, return_counts=True)
+    # print(he)
+    ne = he_u[0][he_u[1]>2]
+
+    return ne
+
diff --git a/test/test_cut_edges.py b/test/test_cut_edges.py
new file mode 100644
index 00000000..a733ac6d
--- /dev/null
+++ b/test/test_cut_edges.py
@@ -0,0 +1,58 @@
+from .context import gpytoolbox as gpy
+from .context import numpy as np
+from .context import unittest
+
+class TestCutEdges(unittest.TestCase):
+
+    def test_two_triangles(self):
+        F = np.array([[0,1,2], [2,1,3]],dtype=int)
+        E = np.array([[1,2]])
+        G,I = gpy.cut_edges(F,E)
+        self.assertTrue(np.all(G==np.array([[0,2,4],[1,3,5]])))
+        self.assertTrue(np.all(I==np.array([0, 2, 1, 1, 2, 3])))
+
+    def test_icosphere(self):
+        V,F = gpy.icosphere(3)
+        E = np.array([[371,573], [573,571], [571,219]])
+        G,I = gpy.cut_edges(F,E)
+
+        b = gpy.boundary_loops(G)
+        print(b)
+        self.assertTrue(b[0].size == 2*E.shape[0])
+        self.assertTrue(np.all(np.sort(b[0]) == np.sort(
+            np.array([
+                557,
+                278,
+                641,
+                437,
+                507,
+                642
+                ]))))
+    
+    def test_bunny(self):
+        V,F = gpy.read_mesh("test/unit_tests_data/bunny_oded.obj")
+        path = np.array([1575,1482,1394,1309,1310,1225,1141])
+        E = np.stack((path[:-1],path[1:]), axis=-1)
+        G,I = gpy.cut_edges(F,E)
+
+        b = gpy.boundary_loops(G)
+        self.assertTrue(b[0].size == 2*E.shape[0])
+        self.assertTrue(np.all(np.sort(b[0]) == np.sort(
+            np.array([
+                 761,
+                 790,
+                 811,
+                 865,
+                 927,
+                 929,
+                 952,
+                1038,
+                2501,
+                2578,
+                2596,
+                2613
+                ]))))
+        
+
+if __name__ == '__main__':
+    unittest.main()
\ No newline at end of file
diff --git a/test/test_non_manifold_edges.py b/test/test_non_manifold_edges.py
new file mode 100644
index 00000000..21e856e0
--- /dev/null
+++ b/test/test_non_manifold_edges.py
@@ -0,0 +1,72 @@
+from .context import gpytoolbox as gpy
+from .context import numpy as np
+from .context import unittest
+
+class TestNonManifoldEdges(unittest.TestCase):
+    # There is not much to test here that goes beyond just inputting the
+    # definition of the function, but we can make sure that a few conditions
+    # are fulfilled.
+
+    def test_single_triangle(self):
+        f = np.array([[0,1,2]],dtype=int)
+        ne = gpy.non_manifold_edges(f)
+        # no non-manifold edges
+        self.assertTrue(len(ne)==0)
+
+    def test_simple_nonmanifold(self):
+        f = np.array([[0,1,2],[0,2,3],[2,0,4]],dtype=int)
+        ne = gpy.non_manifold_edges(f)
+        ne_gt = np.array([[0,2]],dtype=int)
+        self.assertTrue(np.all(ne==ne_gt))
+    
+    def test_bunny(self):
+        _,f = gpy.read_mesh("test/unit_tests_data/bunny_oded.obj")
+        num_faces = f.shape[0]
+        he = gpy.halfedges(f).reshape(-1,2)
+
+        for it in range(100):
+            # pick a random edge in he
+            i = np.random.randint(he.shape[0])
+            random_edge = he[i,:]
+            # now we add a new face that contains this edge
+            new_face = np.array([random_edge[0],random_edge[1],num_faces],dtype=int)
+            # insert this face into a random position in f
+            f_bad = np.insert(f,np.random.randint(f.shape[0]),new_face,axis=0)
+            # are there any non-manifold edges?
+            ne = gpy.non_manifold_edges(f_bad)
+            self.assertTrue(ne.shape[0]==1)
+            # sort random_edge
+            random_edge = np.sort(random_edge)
+            # check that ne is random edge
+            self.assertTrue(np.all(ne==random_edge))
+            # print(random_edge)
+        
+        # now let's add them sequentially
+        f = f_bad.copy()
+        ne_gt = ne.copy()
+        rng = np.random.default_rng(5)
+        for it in range(10):
+            # pick a random edge in he
+            i = rng.integers(he.shape[0])
+            random_edge = he[i,:]
+            # now we add a new face that contains this edge
+            new_face = np.array([random_edge[0],random_edge[1],num_faces+it],dtype=int)
+            # insert this face into a random position in f
+            f = np.insert(f,np.random.randint(f.shape[0]),new_face,axis=0)
+            # are there any non-manifold edges?
+            ne = gpy.non_manifold_edges(f)
+            
+            self.assertTrue(ne.shape[0]==it+2)
+            # sort random_edge
+            random_edge = np.sort(random_edge)
+            ne_gt = np.vstack((ne_gt,random_edge))
+            # sort ne_gt lexicographically
+            ne_gt = np.unique(ne_gt,axis=0,)
+            # should match
+            self.assertTrue(np.all(ne==ne_gt))
+
+        
+
+
+if __name__ == '__main__':
+    unittest.main()
\ No newline at end of file