diff --git a/doc/examples/plot_rag_boundary.py b/doc/examples/plot_rag_boundary.py
new file mode 100644
index 00000000..eeff1589
--- /dev/null
+++ b/doc/examples/plot_rag_boundary.py
@@ -0,0 +1,28 @@
+"""
+==========================
+Region Boundary based RAGs
+==========================
+
+This example demonstrates construction of region boundary based RAGs with the
+`rag_boundary` function.
+"""
+from skimage.future import graph
+from skimage import data, segmentation, color, filters, io
+from matplotlib import pyplot as plt, colors
+
+img = data.coffee()
+gimg = color.rgb2gray(img)
+
+labels = segmentation.slic(img, compactness=30, n_segments=400)
+edges = filters.sobel(gimg)
+edges_rgb = color.gray2rgb(edges)
+
+mimg = segmentation.mark_boundaries(img, labels, (0,0,0))
+
+g = graph.rag_boundary(labels, edges)
+
+cmap = colors.ListedColormap(['#0000ff', '#ff0000'])
+out = graph.draw_rag(labels, g, edges_rgb, node_color="#ffff00", colormap=cmap)
+
+io.imshow(out)
+io.show()
diff --git a/skimage/future/graph/__init__.py b/skimage/future/graph/__init__.py
index 6b7e0aaa..e73c8d19 100644
--- a/skimage/future/graph/__init__.py
+++ b/skimage/future/graph/__init__.py
@@ -1,5 +1,5 @@
 from .graph_cut import cut_threshold, cut_normalized
-from .rag import rag_mean_color, RAG, draw_rag
+from .rag import rag_mean_color, RAG, draw_rag, rag_boundary
 from .graph_merge import merge_hierarchical
 ncut = cut_normalized
 
@@ -9,4 +9,5 @@ __all__ = ['rag_mean_color',
            'ncut',
            'draw_rag',
            'merge_hierarchical',
+           'rag_boundary'
            'RAG']
diff --git a/skimage/future/graph/rag.py b/skimage/future/graph/rag.py
index 6a45447b..6c8fa502 100644
--- a/skimage/future/graph/rag.py
+++ b/skimage/future/graph/rag.py
@@ -311,6 +311,79 @@ def rag_mean_color(image, labels, connectivity=2, mode='distance',
     return graph
 
 
+def rag_boundary(labels, edge_map, connectivity=2):
+    """ Comouter RAG based on region boundaries
+
+    Given an image's initial segmentation and its edge map this method
+    constructs the corresponding Region Adjacency Graph (RAG). Each node in the
+    RAG represents a set of pixels within the image with the same label in
+    `labels`. The weight between two adjacent regions is the average value
+    in `edge_map` along their boundary.
+
+    labels : ndarray
+        The labelled image.
+    edge_map : ndarray
+        This should have the same shape as that of `labels`. For all pixels
+        along the boundary between 2 adjacent regions, the average value of the
+        corresponding pixels in `edge_map` is the edge weight between them.
+    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`.
+
+    Examples
+    --------
+    >>> from skimage import data, segmentation, filters, color
+    >>> from skimage.future import graph
+    >>> img = data.chelsea()
+    >>> labels = segmentation.slic(img)
+    >>> edge_map = filters.sobel(color.rgb2gray(img))
+    >>> rag = graph.rag_mean_color(labels, edge_map)
+
+    """
+
+    graph = RAG()
+
+    #Computing the relative indices of the neighbors
+    nbr_indices = list(np.ndindex(*[2]*labels.ndim))
+    del nbr_indices[0]
+    nbr_indices_arr = ([idx for idx in nbr_indices if np.linalg.norm(idx)
+                   <= connectivity])
+
+
+    iter_shape = tuple(np.array(labels.shape) - 1)
+
+    for index in np.ndindex(iter_shape):
+
+        index_arr = np.array(index)
+        current = labels[index]
+        graph.add_node(current, {'labels':[current]})
+
+        for nbr_index in nbr_indices_arr:
+
+            adjacent_idx = tuple(index_arr + nbr_index)
+            adjacent = labels[adjacent_idx]
+
+            if current==adjacent:
+                continue
+
+            if graph.has_edge(current, adjacent):
+                graph[current][adjacent]['pixel count'] += 2
+                intensity = edge_map[index] + edge_map[adjacent_idx]
+                graph[current][adjacent]['total intensity'] += intensity
+            else:
+                graph.add_edge(current, adjacent)
+                graph[current][adjacent]['pixel count'] = 2
+                intensity = edge_map[index] + edge_map[adjacent_idx]
+                graph[current][adjacent]['total intensity'] = intensity
+
+    for (x, y, data) in graph.edges_iter(data=True):
+        data['weight'] = data['total intensity']/data['pixel count']
+
+    return graph
+
+
 def draw_rag(labels, rag, img, border_color=None, node_color='#ffff00',
              edge_color='#00ff00', colormap=None, thresh=np.inf,
              desaturate=False, in_place=True):
diff --git a/skimage/future/graph/tests/test_rag.py b/skimage/future/graph/tests/test_rag.py
index e5882e20..8f1df909 100644
--- a/skimage/future/graph/tests/test_rag.py
+++ b/skimage/future/graph/tests/test_rag.py
@@ -196,3 +196,20 @@ def test_generic_rag_3d():
     assert h.has_edge(0, 1) and h.has_edge(0, 3) and not h.has_edge(0, 7)
     k = graph.RAG(labels, connectivity=3)
     assert k.has_edge(0, 1) and k.has_edge(1, 2) and k.has_edge(2, 5)
+
+
+def test_rag_boundary():
+    labels = np.zeros((16, 16), dtype='uint8')
+    edge_map = np.zeros_like(labels, dtype=float)
+
+    edge_map[8,:] = 0.5
+    edge_map[:,8] = 1.0
+
+    labels[:8, :8] = 1
+    labels[:8, 8:] = 2
+    labels[8:, :8] = 3
+    labels[8:, 8:] = 4
+
+    g = graph.rag_boundary(labels, edge_map, connectivity=1)
+    assert len(g.nodes()) == 4
+    assert len(g.edges()) == 4
\ No newline at end of file