diff --git a/.travis.yml b/.travis.yml
index 169b512b..1f292dfa 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -54,6 +54,7 @@ before_install:
     - pip install cython
     - pip install flake8
     - pip install six
+    - pip install networkx
 
     - pip install nose-cov
     - pip install coveralls
diff --git a/doc/examples/plot_rag.py b/doc/examples/plot_rag.py
new file mode 100644
index 00000000..b26a266d
--- /dev/null
+++ b/doc/examples/plot_rag.py
@@ -0,0 +1,81 @@
+"""
+=======================
+Region Adjacency Graphs
+=======================
+
+This example demonstrates the use of the `merge_nodes` function of a Region
+Adjacency Graph (RAG). The `RAG` class represents a undirected weighted graph
+which inherits from `networkx.graph` class. When a new node is formed by
+merging two nodes, the edge weight of all the edges incident on the resulting
+node can be updated by a user defined function `weight_func`.
+
+The default behaviour is to use the smaller edge weight in case of a conflict.
+The example below also shows how to use a custom function to select the larger
+weight instead.
+
+"""
+from skimage.graph import rag
+import networkx as nx
+from matplotlib import pyplot as plt
+import numpy as np
+
+
+def max_edge(g, src, dst, n):
+    """Callback to handle merging nodes by choosing maximum weight.
+
+    Returns either the weight between (`src`, `n`) or (`dst`, `n`)
+    in `g` or the maximum of the two when both exist.
+
+    Parameters
+    ----------
+    g : RAG
+        The graph under consideration.
+    src, dst : int
+        The verices in `g` to be merged.
+    n : int
+        A neighbor of `src` or `dst` or both.
+
+    Returns
+    -------
+    weight : float
+        The weight between (`src`, `n`) or (`dst`, `n`) in `g` or the
+        maximum of the two when both exist.
+
+    """
+
+    w1 = g[n].get(src, {'weight': -np.inf})['weight']
+    w2 = g[n].get(dst, {'weight': -np.inf})['weight']
+    return max(w1, w2)
+
+
+def display(g, title):
+    """Displays a graph with the given title."""
+    pos = nx.circular_layout(g)
+    plt.figure()
+    plt.title(title)
+    nx.draw(g, pos)
+    nx.draw_networkx_edge_labels(g, pos, font_size=20)
+
+
+g = rag.RAG()
+g.add_edge(1, 2, weight=10)
+g.add_edge(2, 3, weight=20)
+g.add_edge(3, 4, weight=30)
+g.add_edge(4, 1, weight=40)
+g.add_edge(1, 3, weight=50)
+
+# Assigning dummy labels.
+for n in g.nodes():
+    g.node[n]['labels'] = [n]
+
+gc = g.copy()
+
+display(g, "Original Graph")
+
+g.merge_nodes(1, 3)
+display(g, "Merged with default (min)")
+
+gc.merge_nodes(1, 3, weight_func=max_edge)
+display(gc, "Merged with max")
+
+plt.show()
diff --git a/doc/examples/plot_rag_mean_color.py b/doc/examples/plot_rag_mean_color.py
new file mode 100644
index 00000000..7a1d24b9
--- /dev/null
+++ b/doc/examples/plot_rag_mean_color.py
@@ -0,0 +1,29 @@
+"""
+================
+RAG Thresholding
+================
+
+This example constructs a Region Adjacency Graph (RAG) and merges regions
+which are similar in color. We construct a RAG and define edges as the
+difference in mean color. We then join regions with similar mean color.
+
+"""
+
+from skimage import graph, data, io, segmentation, color
+from matplotlib import pyplot as plt
+
+
+img = data.coffee()
+
+labels1 = segmentation.slic(img, compactness=30, n_segments=400)
+out1 = color.label2rgb(labels1, img, kind='avg')
+
+g = graph.rag_mean_color(img, labels1)
+labels2 = graph.cut_threshold(labels1, g, 30)
+out2 = color.label2rgb(labels2, img, kind='avg')
+
+plt.figure()
+io.imshow(out1)
+plt.figure()
+io.imshow(out2)
+io.show()
diff --git a/requirements.txt b/requirements.txt
index 48e03b8b..7543d867 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -2,3 +2,4 @@ cython>=0.17
 matplotlib>=1.0
 numpy>=1.6
 six>=1.3.0
+networkx>=1.8.0
diff --git a/skimage/graph/__init__.py b/skimage/graph/__init__.py
index a335971d..68c1206c 100644
--- a/skimage/graph/__init__.py
+++ b/skimage/graph/__init__.py
@@ -1,9 +1,14 @@
 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
 
 __all__ = ['shortest_path',
            'MCP',
            'MCP_Geometric',
            'MCP_Connect',
            'MCP_Flexible',
-           'route_through_array']
\ No newline at end of file
+           'route_through_array',
+           'rag_mean_color',
+           'cut_threshold',
+           'RAG']
diff --git a/skimage/graph/graph_cut.py b/skimage/graph/graph_cut.py
new file mode 100644
index 00000000..91d64aec
--- /dev/null
+++ b/skimage/graph/graph_cut.py
@@ -0,0 +1,58 @@
+import networkx as nx
+import numpy as np
+
+
+def cut_threshold(labels, rag, thresh):
+    """Combine regions seperated by weight less than threshold.
+
+    Given an image's labels and its RAG, output new labels by
+    combining regions whose nodes are seperated by a weight less
+    than the given threshold.
+
+    Parameters
+    ----------
+    labels : ndarray
+        The array of labels.
+    rag : RAG
+        The region adjacency graph.
+    thresh : float
+        The threshold. Regions connected by edges with smaller weights are
+        combined.
+
+    Returns
+    -------
+    out : ndarray
+        The new labelled array.
+
+    Examples
+    --------
+    >>> from skimage import data, graph, segmentation
+    >>> img = data.lena()
+    >>> labels = segmentation.slic(img)
+    >>> rag = graph.rag_mean_color(img, labels)
+    >>> new_labels = graph.cut_threshold(labels, rag, 10)
+
+    References
+    ----------
+    .. [1] Alain Tremeau and Philippe Colantoni
+           "Regions Adjacency Graph Applied To Color Image Segmentation"
+           http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.11.5274
+
+    """
+    # Because deleting edges while iterating through them produces an error.
+    to_remove = [(x, y) for x, y, d in rag.edges_iter(data=True)
+                 if d['weight'] >= thresh]
+    rag.remove_edges_from(to_remove)
+
+    comps = nx.connected_components(rag)
+
+    # We construct an array which can map old labels to the new ones.
+    # All the labels within a connected component are assigned to a single
+    # label in the output.
+    map_array = np.arange(labels.max() + 1, dtype=labels.dtype)
+    for i, nodes in enumerate(comps):
+        for node in nodes:
+            for label in rag.node[node]['labels']:
+                map_array[label] = i
+
+    return map_array[labels]
diff --git a/skimage/graph/rag.py b/skimage/graph/rag.py
new file mode 100644
index 00000000..84721736
--- /dev/null
+++ b/skimage/graph/rag.py
@@ -0,0 +1,203 @@
+import networkx as nx
+import numpy as np
+from scipy.ndimage import filters
+from scipy import ndimage as nd
+
+
+def min_weight(graph, src, dst, n):
+    """Callback to handle merging nodes by choosing minimum weight.
+
+    Returns either the weight between (`src`, `n`) or (`dst`, `n`)
+    in `graph` or the minumum of the two when both exist.
+
+    Parameters
+    ----------
+    graph : RAG
+        The graph under consideration.
+    src, dst : int
+        The verices in `graph` to be merged.
+    n : int
+        A neighbor of `src` or `dst` or both.
+
+    Returns
+    -------
+    weight : float
+        The weight between (`src`, `n`) or (`dst`, `n`) in `graph` or the
+        minumum of the two when both exist.
+
+    """
+
+    # cover the cases where n only has edge to either `src` or `dst`
+    default = {'weight': np.inf}
+    w1 = graph[n].get(src, default)['weight']
+    w2 = graph[n].get(dst, default)['weight']
+    return min(w1, w2)
+
+
+class RAG(nx.Graph):
+
+    """
+    The Region Adjacency Graph (RAG) of an image, subclasses
+    `networx.Graph <http://networkx.github.io/documentation/latest/reference/classes.graph.html>`_
+    """
+
+    def merge_nodes(self, src, dst, weight_func=min_weight, extra_arguments=[],
+                    extra_keywords={}):
+        """Merge node `src` into `dst`.
+
+        The new combined node is adjacent to all the neighbors of `src`
+        and `dst`. `weight_func` is called to decide the weight of edges
+        incident on the new node.
+
+        Parameters
+        ----------
+        src, dst : int
+            Nodes to be merged.
+        weight_func : callable, optional
+            Function to decide edge weight of edges incident on the new node.
+            For each neighbor `n` for `src and `dst`, `weight_func` will be
+            called as follows: `weight_func(src, dst, n, *extra_arguments,
+            **extra_keywords)`. `src`, `dst` and `n` are IDs of vertices in the
+            RAG object which is in turn a subclass of
+            `networkx.Graph`.
+        extra_arguments : sequence, optional
+            The sequence of extra positional arguments passed to
+            `weight_func`.
+        extra_keywords : dictionary, optional
+            The dict of keyword arguments passed to the `weight_func`.
+
+        """
+        src_nbrs = set(self.neighbors(src))
+        dst_nbrs = set(self.neighbors(dst))
+        neighbors = (src_nbrs & dst_nbrs) - set([src, dst])
+
+        for neighbor in neighbors:
+            w = weight_func(self, src, dst, neighbor, *extra_arguments,
+                            **extra_keywords)
+            self.add_edge(neighbor, dst, weight=w)
+
+        self.node[dst]['labels'] += self.node[src]['labels']
+        self.remove_node(src)
+
+
+def _add_edge_filter(values, graph):
+    """Create edge in `g` between the first element of `values` and the rest.
+
+    Add an edge between the first element in `values` and
+    all other elements of `values` in the graph `g`. `values[0]`
+    is expected to be the central value of the footprint used.
+
+    Parameters
+    ----------
+    values : array
+        The array to process.
+    graph : RAG
+        The graph to add edges in.
+
+    Returns
+    -------
+    0 : int
+        Always returns 0. The return value is required so that `generic_filter`
+        can put it in the output array.
+
+    """
+    values = values.astype(int)
+    current = values[0]
+    for value in values[1:]:
+        graph.add_edge(current, value)
+
+    return 0
+
+
+def rag_mean_color(image, labels, connectivity=2):
+    """Compute the Region Adjacency Graph using mean colors.
+
+    Given an image and its initial segmentation, this method constructs the
+    corresponsing Region Adjacency Graph (RAG). Each node in the RAG
+    represents a set of pixels within `image` with the same label in `labels`.
+    The weight between two adjacent regions is the difference in their mean
+    color.
+
+    Parameters
+    ----------
+    image : ndarray, shape(M, N, [..., P,] 3)
+        Input image.
+    labels : ndarray, shape(M, N, [..., P,])
+        The labelled image. This should have one dimension less than
+        `image`. If `image` has dimensions `(M, N, 3)` `labels` should have
+         dimensions `(M, N)`.
+    connectivity : int, optional
+        Pixels with a squared distance less than `connectivity` from each other
+        are considered adjacent. It can range from 1 to `labels.ndim`. Its
+        behavior is the same as `connectivity` parameter in
+        `scipy.ndimage.filters.generate_binary_structure`.
+
+    Returns
+    -------
+    out : RAG
+        The region adjacency graph.
+
+    Examples
+    --------
+    >>> from skimage import data, graph, segmentation
+    >>> img = data.lena()
+    >>> labels = segmentation.slic(img)
+    >>> rag = graph.rag_mean_color(img, labels)
+
+    References
+    ----------
+    .. [1] Alain Tremeau and Philippe Colantoni
+           "Regions Adjacency Graph Applied To Color Image Segmentation"
+           http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.11.5274
+
+    """
+    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)
+
+    # For example
+    # if labels.ndim = 2 and connectivity = 1
+    # fp = [[0,0,0],
+    #       [0,1,1],
+    #       [0,1,0]]
+    #
+    # if labels.ndim = 2 and connectivity = 2
+    # fp = [[0,0,0],
+    #       [0,1,1],
+    #       [0,1,1]]
+
+    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']
+        graph[x][y]['weight'] = np.linalg.norm(diff)
+
+    return graph
diff --git a/skimage/graph/tests/test_rag.py b/skimage/graph/tests/test_rag.py
new file mode 100644
index 00000000..af7481eb
--- /dev/null
+++ b/skimage/graph/tests/test_rag.py
@@ -0,0 +1,64 @@
+import numpy as np
+from skimage import graph
+
+
+def max_edge(g, src, dst, n):
+    default = {'weight': -np.inf}
+    w1 = g[n].get(src, default)['weight']
+    w2 = g[n].get(dst, default)['weight']
+    return max(w1, w2)
+
+
+def test_rag_merge():
+    g = graph.rag.RAG()
+
+    for i in range(5):
+        g.add_node(i, {'labels': [i]})
+
+    g.add_edge(0, 1, {'weight': 10})
+    g.add_edge(1, 2, {'weight': 20})
+    g.add_edge(2, 3, {'weight': 30})
+    g.add_edge(3, 0, {'weight': 40})
+    g.add_edge(0, 2, {'weight': 50})
+    g.add_edge(3, 4, {'weight': 60})
+
+    gc = g.copy()
+
+    # We merge nodes and ensure that the minimum weight is chosen
+    # when there is a conflict.
+    g.merge_nodes(0, 2)
+    assert g.edge[1][2]['weight'] == 10
+    assert g.edge[2][3]['weight'] == 30
+
+    # We specify `max_edge` as `weight_func` as ensure that maximum
+    # weight is chosen in case on conflict
+    gc.merge_nodes(0, 2, weight_func=max_edge)
+    assert gc.edge[1][2]['weight'] == 20
+    assert gc.edge[2][3]['weight'] == 40
+
+    g.merge_nodes(1, 4)
+    g.merge_nodes(2, 3)
+    g.merge_nodes(3, 4)
+    assert sorted(g.node[4]['labels']) == list(range(5))
+    assert g.edges() == []
+
+
+def test_threshold_cut():
+
+    img = np.zeros((100, 100, 3), dtype='uint8')
+    img[:50, :50] = 255, 255, 255
+    img[:50, 50:] = 254, 254, 254
+    img[50:, :50] = 2, 2, 2
+    img[50:, 50:] = 1, 1, 1
+
+    labels = np.zeros((100, 100), dtype='uint8')
+    labels[:50, :50] = 0
+    labels[:50, 50:] = 1
+    labels[50:, :50] = 2
+    labels[50:, 50:] = 3
+
+    rag = graph.rag_mean_color(img, labels)
+    new_labels = graph.cut_threshold(labels, rag, 10)
+
+    # Two labels
+    assert new_labels.max() == 1