diff --git a/skimage/graph/__init__.py b/skimage/graph/__init__.py
index 2d774b44..aca73bc4 100644
--- a/skimage/graph/__init__.py
+++ b/skimage/graph/__init__.py
@@ -2,6 +2,7 @@ from .spath import shortest_path
 from .mcp import MCP, MCP_Geometric, MCP_Connect, MCP_Flexible, route_through_array
 from .rag import rag_mean_color, RAG
 from .graph_cut import cut_threshold, cut_normalized
+ncut = cut_normalized
 
 __all__ = ['shortest_path',
            'MCP',
@@ -12,4 +13,5 @@ __all__ = ['shortest_path',
            'rag_mean_color',
            'cut_threshold',
            'cut_normalized',
+           'ncut',
            'RAG']
diff --git a/skimage/graph/_ncut.py b/skimage/graph/_ncut.py
index 4874b252..fbefca27 100644
--- a/skimage/graph/_ncut.py
+++ b/skimage/graph/_ncut.py
@@ -21,20 +21,20 @@ def DW_matrices(graph):
         The weight matrix of the graph. `W[i, j]` is the weight of the edge
         joining `i` to `j`.
     """
-    #Cause sparse.eigsh prefers CSC format
+    # sparse.eighsh is most efficient with CSC-formatted input
     W = nx.to_scipy_sparse_matrix(graph, format='csc')
     entries = W.sum(axis=0)
     D = sparse.dia_matrix((entries, 0), shape=W.shape).tocsc()
     return D, W
 
 
-def ncut_cost(mask, D, W):
+def ncut_cost(cut, D, W):
     """Returns the N-cut cost of a bi-partition of a graph.
 
     Parameters
     ----------
-    mask : ndarray
-        The mask for the nodes in the graph. Nodes corrsesponding to a `True`
+    cut : ndarray
+        The mask for the nodes in the graph. Nodes corressponding to a `True`
         value are in one set.
     D : csc_matrix
         The diagonal matrix of the graph.
@@ -45,15 +45,22 @@ def ncut_cost(mask, D, W):
     -------
     cost : float
         The cost of performing the N-cut.
+
+    References
+    ----------
+    .. [1] Normalized Cuts and Image Segmentation, Jianbo Shi and
+           Jitendra Malik, IEEE Transactions on Pattern Analysis and Machine
+           Intelligence, Page 889, Equation 2.
     """
-    mask = np.array(mask)
-    cut = _ncut_cy.cut_cost(mask, W)
+    cut = np.array(cut)
+    cut_cost = _ncut_cy.cut_cost(cut, W)
 
-    # Cause D has elements only along diagonal
-    assoc_a = D.data[mask].sum()
-    assoc_b = D.data[np.logical_not(mask)].sum()
+    # D has elements only along the diagonal, one per node, so we can directly
+    # index the data attribute with cut.
+    assoc_a = D.data[cut].sum()
+    assoc_b = D.data[np.logical_not(cut)].sum()
 
-    return (cut / assoc_a) + (cut / assoc_b)
+    return (cut_cost / assoc_a) + (cut_cost / assoc_b)
 
 
 def normalize(a):
diff --git a/skimage/graph/_ncut_cy.pyx b/skimage/graph/_ncut_cy.pyx
index 8c416cbc..0f859edf 100644
--- a/skimage/graph/_ncut_cy.pyx
+++ b/skimage/graph/_ncut_cy.pyx
@@ -16,52 +16,65 @@ def argmin2(cnp.float64_t[:] array):
 
     Returns
     -------
-    i : int
+    min_idx2 : int
         The index of the second smallest value.
     """
     cdef cnp.float64_t min1 = np.inf
     cdef cnp.float64_t min2 = np.inf
-    cdef Py_ssize_t i1 = 0
-    cdef Py_ssize_t i2 = 0
+    cdef Py_ssize_t min_idx1 = 0
+    cdef Py_ssize_t min_idx2 = 0
     cdef Py_ssize_t i = 0
+    cdef Py_ssize_t n = array.shape[0]
 
-    while i < array.shape[0]:
+    for i in range(n):
         x = array[i]
         if x < min1:
             min2 = min1
-            i2 = i1
+            min_idx2 = min_idx1
             min1 = x
-            i1 = i
+            min_idx1 = i
         elif x > min1 and x < min2:
             min2 = x
-            i2 = i
+            min_idx2 = i
         i += 1
 
-    return i2
+    return min_idx2
 
 
-def cut_cost(mask, W):
-    mask = np.array(mask)
+def cut_cost(cut, W):
+    """Return the total weight of crossing edges in a bi-partition.
 
+    Parameters
+    ----------
+    cut : array
+        A array of booleans. Elements set to `True` belong to one
+        set.
+    W : array
+        The weight matrix of the graph.
+
+    Returns
+    -------
+    cost : float
+        The total weight of crossing edges.
+    """
+    cdef cnp.ndarray[cnp.uint8_t, cast = True] cut_mask = np.array(cut)
     cdef Py_ssize_t num_rows, num_cols
     cdef cnp.int32_t row, col
     cdef cnp.int32_t[:] indices = W.indices
     cdef cnp.int32_t[:] indptr = W.indptr
-    cdef cnp.float64_t[:] data = W.data
+    cdef cnp.double_t[:] data = W.data.astype(np.double)
     cdef cnp.int32_t row_index
     cdef cnp.double_t cost = 0
 
     num_rows = W.shape[0]
     num_cols = W.shape[1]
 
-    col = 0
-    while col < num_cols:
-        row_index = indptr[col]
-        while row_index < indptr[col+1]:
+    for col in range(num_cols):
+        for row_index in range(indptr[col], indptr[col + 1]):
             row = indices[row_index]
-            if mask[row] != mask[col]:
-                cost += <cnp.double_t>data[row_index]
+            if cut_mask[row] != cut_mask[col]:
+                cost += data[row_index]
             row_index += 1
         col += 1
 
-    return cost*0.5
+    return cost * 0.5
diff --git a/skimage/graph/graph_cut.py b/skimage/graph/graph_cut.py
index af5048fd..9c133bd3 100644
--- a/skimage/graph/graph_cut.py
+++ b/skimage/graph/graph_cut.py
@@ -111,6 +111,65 @@ def cut_normalized(labels, rag, thresh=0.001, num_cuts=10):
     return map_array[labels]
 
 
+def partition_by_cut(cut, rag):
+    """Compute resulting subgraphs from given bi-parition.
+
+    Parameters
+    ----------
+    cut : array
+        A array of booleans. Elements set to `True` belong to one
+        set.
+    rag : RAG
+        The Region Adjacency Graph.
+
+    Returns
+    -------
+    sub1, sub2 : RAG
+        The two resulting subgraphs from the bi-partition.
+    """
+    nodes1 = [n for i, n in enumerate(rag.nodes()) if cut[i]]
+    nodes2 = [n for i, n in enumerate(rag.nodes()) if not cut[i]]
+
+    sub1 = rag.subgraph(nodes1)
+    sub2 = rag.subgraph(nodes2)
+
+    return sub1, sub2
+
+
+def get_min_ncut(ev, d, w, num_cuts):
+    """Threshold an eigen vector evenly, to determine minimum ncut.
+
+    Parameters
+    ----------
+    ev : array
+        The eigenvector to threshold.
+    d : ndarray
+        The diagonal matrix of the graph.
+    w : ndarray
+        The weight matrix of the graph.
+    num_cuts : int
+        The number of evenly spaced thresholds to check for.
+
+    Returns
+    -------
+    threshold, mcut : float
+        The threshold which produced the minimum ncut, and the value of the
+        ncut itself.
+    """
+    mcut = np.inf
+
+    # Refer Shi & Malik 2001, Section 3.1.3, Page 892
+    # Perform evenly spaced n-cuts and determine the optimal one.
+    for t in np.linspace(0, 1, num_cuts, endpoint=False):
+        mask = ev > t
+        cost = _ncut.ncut_cost(mask, d, w)
+        if cost < mcut:
+            mcut = cost
+            threshold = t
+
+    return threshold, mcut
+
+
 def _ncut_relabel(rag, thresh, num_cuts, map_array):
     """Perform Normalized Graph cut on the Region Adjacency Graph.
 
@@ -135,58 +194,41 @@ def _ncut_relabel(rag, thresh, num_cuts, map_array):
         the function.
     """
     d, w = _ncut.DW_matrices(rag)
-    error = False
+    stop = False
+    m = w.shape[0]
 
-    try:
-        m = w.shape[0]
+    if m > 2:
         d2 = d.copy()
         # Since d is diagonal, we can directly operate on it's data
         # the inverse
-        d2.data = 1.0/d2.data
+        d2.data = 1.0 / d2.data
         # the square root
         d2.data = np.sqrt(d2.data)
-        # Refer to Equation 7
-        vals, vectors = linalg.eigsh(d2*(d - w)*d2, which='SM',
+        # Refer Shi & Malik 2001, Equation 7, Page 891
+        vals, vectors = linalg.eigsh(d2 * (d - w) * d2, which='SM',
                                      k=min(100, m - 2))
-    except ValueError:
-        # k is too less, happens when the graph is of size 1
-        error = True
+    else:
+        stop = True
 
-    if not error:
-        # Refer Section 3.2.3
+    if not stop:
+        # Pick second smalles eigen vector.
+        # Refer Shi & Malik 2001, Section 3.2.3, Page 893
         vals, vectors = np.real(vals), np.real(vectors)
         index2 = _ncut_cy.argmin2(vals)
+        ev = _ncut.normalize(vectors[:, index2])
 
-        ev = np.real(vectors[:, index2])
-        ev = _ncut.normalize(ev)
-
-        mcut = np.inf
-        threshold = None
-        # Refer Section 3.1.3
-        # Perform evenly spaced n-cuts and determine the optimal one.
-        for t in np.linspace(0, 1, num_cuts, endpoint=False):
-            mask = ev > t
-            cost = _ncut.ncut_cost(mask, d, w)
-            if cost < mcut:
-                mcut = cost
-                threshold = t
-
+        threshold, mcut = get_min_ncut(ev, d, w, num_cuts)
         if (mcut < thresh):
-            mask = ev > threshold
-
-            nodes1 = [n for i, n in enumerate(rag.nodes()) if mask[i]]
-            nodes2 = [n for i, n in enumerate(rag.nodes()) if not mask[i]]
-
+            cut_mask = ev > threshold
             # Sub divide and perform N-cut again
-            sub1 = rag.subgraph(nodes1)
-            sub2 = rag.subgraph(nodes2)
+            # Refer Shi & Malik 2001, Section 3.2.5, Page 893
+            sub1, sub2 = partition_by_cut(cut_mask, rag)
 
-            # Refer Section 3.2.5
             _ncut_relabel(sub1, thresh, num_cuts, map_array)
             _ncut_relabel(sub2, thresh, num_cuts, map_array)
             return
 
-    # Either an errornous condition occurred, or N-cut wasn't small enough.
+    # The N-cut wasn't small enough, or could not be computed.
     # The remaining graph is a region.
     # Assign `ncut label` by picking any label from the existing nodes, since
     # `labels` are unique, 'ncut label' is also unique.