diff --git a/skimage/graph/__init__.py b/skimage/graph/__init__.py
index 68c1206c..9fe00d3c 100644
--- a/skimage/graph/__init__.py
+++ b/skimage/graph/__init__.py
@@ -1,7 +1,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
+from .graph_cut import cut_threshold, cut_n
 
 __all__ = ['shortest_path',
            'MCP',
@@ -11,4 +11,5 @@ __all__ = ['shortest_path',
            'route_through_array',
            'rag_mean_color',
            'cut_threshold',
+           'cut_n',
            'RAG']
diff --git a/skimage/graph/_ncut.py b/skimage/graph/_ncut.py
new file mode 100644
index 00000000..9b321bd8
--- /dev/null
+++ b/skimage/graph/_ncut.py
@@ -0,0 +1,38 @@
+import networkx as nx
+import numpy as np
+from scipy import sparse
+
+def DW_matrix(graph):
+    #print graph[4][4]['weight']
+    W = nx.to_scipy_sparse_matrix(graph, format='csc')
+    entries = W.sum(0)
+    D = sparse.dia_matrix( (entries,0),shape = W.shape).tocsc()
+    #print W[4,4]
+    
+    m,n = W.shape
+    #for i in range(n):
+    #    W[i,i] = 1.0
+    return D,W
+
+def ncut_cost(mask,D,W):
+
+    mask = np.array(mask)
+    mask_list = [ np.logical_xor(mask[i], mask) for i in range(mask.shape[0])]
+    mask_array = np.array(mask_list)
+       
+    cut = float(W[mask_array].sum()/2.0)
+    #print W.todense()
+    #print mask_array.astype(int)
+    #print "cut=",cut
+    
+    assoc_a = D.data[mask].sum()
+    assoc_b = D.data[np.logical_not(mask)].sum()
+    
+    #print cut
+    #print assoc_a,assoc_b
+    return (cut/assoc_a) + (cut/assoc_b)
+
+def norml(a):
+    mi = a.min()
+    mx = a.max()
+    return (a-mi)/(mx-mi)
diff --git a/skimage/graph/_ncut_cy.pyx b/skimage/graph/_ncut_cy.pyx
new file mode 100644
index 00000000..fc3c0f3b
--- /dev/null
+++ b/skimage/graph/_ncut_cy.pyx
@@ -0,0 +1,27 @@
+# cython: cdivision=True
+# cython: boundscheck=False
+# cython: nonecheck=False
+# cython: wraparound=False
+cimport numpy as cnp
+import numpy as np
+
+def argmin2(cnp.float64_t[:] array):
+    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 i = 0
+
+    while i < array.shape[0]:
+        x = array[i]
+        if x < min1 :
+            min2 = min1
+            i2 = i1
+            min1 = x
+            i1 = i
+        elif x > min1 and x < min2 :
+            min2 = x
+            i2 = i
+        i += 1
+        
+    return i2
diff --git a/skimage/graph/graph_cut.py b/skimage/graph/graph_cut.py
index d9fcfdef..15f6ff34 100644
--- a/skimage/graph/graph_cut.py
+++ b/skimage/graph/graph_cut.py
@@ -4,7 +4,11 @@ except ImportError:
     import warnings
     warnings.warn('"cut_threshold" requires networkx')
 import numpy as np
-
+import _ncut
+import _ncut_cy
+from scipy.sparse import linalg
+from scipy.sparse.linalg.eigen.arpack.arpack import ArpackNoConvergence as ANC
+from scipy.sparse.linalg.eigen.arpack.arpack import ArpackError as APE
 
 def cut_threshold(labels, rag, thresh):
     """Combine regions seperated by weight less than threshold.
@@ -60,3 +64,67 @@ def cut_threshold(labels, rag, thresh):
                 map_array[label] = i
 
     return map_array[labels]
+
+
+def cut_n(labels, rag, thresh):
+
+    _ncut_relabel(rag,thresh)
+    
+    from_ = range(labels.max()+1)
+    to = [ rag.node[x]['ncut label']  for x in from_ ]
+    map_array = np.array(to)
+    
+    return map_array[labels]
+
+def _ncut_relabel(rag, cut_thresh = 0.0001):
+    d, w = _ncut.DW_matrix(rag)
+    error = False
+
+    try:
+        m = w.shape[0]
+        vals,vectors = linalg.eigsh(d-w,M=d,which='SM',k = min(100,m-2))
+    except ANC as e:
+        vals = e.eigenvalues
+        vectors = e.eigenvectors
+        if len(vals) == 0:
+            error = True
+    except ValueError:
+        error = True
+    except APE:
+        error = True
+    
+    if not error :
+        vals,vectors = np.real(vals), np.real(vectors)
+        index2 = _ncut_cy.argmin2(vals)
+
+        ev = np.real(vectors[:,index2])
+        ev = _ncut.norml(ev)
+
+        mcut = np.inf
+        thresh = None
+        for t in np.arange(0,1,0.1):
+            mask = ev > t
+            cost = _ncut.ncut_cost(mask,d,w)
+            if cost < mcut :
+                mcut = cost
+                thresh = t
+
+        if ( mcut < cut_thresh ):
+            mask = ev > thresh
+
+            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]]
+
+            sub1 = rag.subgraph(nodes1)
+            sub2 = rag.subgraph(nodes2)
+
+            _ncut_relabel(sub1,cut_thresh)
+            _ncut_relabel(sub2, cut_thresh)
+            return
+
+    node = rag.nodes()[0]
+    new_label = rag.node[node]['labels'][0]
+    for n in rag.nodes():
+        rag.node[n]['ncut label'] = new_label
+
+    
diff --git a/skimage/graph/rag.py b/skimage/graph/rag.py
index 7c1b7f07..4592ea36 100644
--- a/skimage/graph/rag.py
+++ b/skimage/graph/rag.py
@@ -11,6 +11,7 @@ except ImportError:
 import numpy as np
 from scipy.ndimage import filters
 from scipy import ndimage as nd
+import math
 
 
 def min_weight(graph, src, dst, n):
@@ -210,3 +211,51 @@ def rag_mean_color(image, labels, connectivity=2):
         graph[x][y]['weight'] = np.linalg.norm(diff)
 
     return graph
+ 
+def rag_similarity(image, labels, connectivity=2, sigma = 30.0):
+    graph = RAG()
+
+    # The footprint is constructed in such a way that the first
+    # element in the array being passed to _add_edge_filter is
+    # the central value.
+    fp = nd.generate_binary_structure(labels.ndim, connectivity)
+    for d in range(fp.ndim):
+        fp = fp.swapaxes(0, d)
+        fp[0, ...] = 0
+        fp = fp.swapaxes(0, d)
+
+
+    filters.generic_filter(
+        labels,
+        function=_add_edge_filter,
+        footprint=fp,
+        mode='nearest',
+        output=np.zeros(labels.shape, dtype=np.uint8),
+        extra_arguments=(graph,))
+
+    for n in graph:
+        graph.node[n].update({'labels': [n],
+                              'pixel count': 0,
+                              'total color': np.array([0, 0, 0],
+                                                      dtype=np.double)})
+
+    for index in np.ndindex(labels.shape):
+        current = labels[index]
+        graph.node[current]['pixel count'] += 1
+        graph.node[current]['total color'] += image[index]
+
+    for n in graph:
+        graph.node[n]['mean color'] = (graph.node[n]['total color'] /
+                                       graph.node[n]['pixel count'])
+
+    for x, y in graph.edges_iter():
+        diff = graph.node[x]['mean color'] - graph.node[y]['mean color']
+        diff = np.linalg.norm(diff)
+    
+        #if diff == 0:
+        #    graph[x][y]['weight'] = 99999
+        #else:
+        graph[x][y]['weight'] = math.e**(-(diff**2)/sigma)
+        #print graph[x][y]['weight']
+        #print "diff",diff,"weight",math.e**(-(diff**2)/sigma)
+    return graph
diff --git a/skimage/graph/setup.py b/skimage/graph/setup.py
index 463d2739..edc7b653 100644
--- a/skimage/graph/setup.py
+++ b/skimage/graph/setup.py
@@ -17,6 +17,7 @@ def configuration(parent_package='', top_path=None):
     cython(['_spath.pyx'], working_path=base_path)
     cython(['_mcp.pyx'], working_path=base_path)
     cython(['heap.pyx'], working_path=base_path)
+    cython(['_ncut_cy.pyx'], working_path=base_path)
 
     config.add_extension('_spath', sources=['_spath.c'],
                          include_dirs=[get_numpy_include_dirs()])
@@ -24,7 +25,8 @@ def configuration(parent_package='', top_path=None):
                          include_dirs=[get_numpy_include_dirs()])
     config.add_extension('heap', sources=['heap.c'],
                          include_dirs=[get_numpy_include_dirs()])
-
+    config.add_extension('_ncut_cy', sources=['_ncut_cy.c'],
+                         include_dirs=[get_numpy_include_dirs()])
     return config
 
 if __name__ == '__main__':