diff --git a/doc/examples/plot_join_segmentations.py b/doc/examples/plot_join_segmentations.py
new file mode 100644
index 00000000..02600ba6
--- /dev/null
+++ b/doc/examples/plot_join_segmentations.py
@@ -0,0 +1,69 @@
+"""
+==========================================
+Find the intersection of two segmentations
+==========================================
+
+When segmenting an image, you may want to combine multiple alternative
+segmentations. The `skimage.segmentation.join_segmentations` function
+computes the join of two segmentations, in which a pixel is placed in
+the same segment if and only if it is in the same segment in _both_
+segmentations.
+"""
+
+import numpy as np
+from scipy import ndimage as nd
+import matplotlib.pyplot as plt
+import matplotlib as mpl
+
+from skimage.filter import sobel
+from skimage.segmentation import slic, join_segmentations
+from skimage.morphology import watershed
+
+from skimage import data
+
+coins = data.coins()
+
+# make segmentation using edge-detection and watershed
+edges = sobel(coins)
+markers = np.zeros_like(coins)
+foreground, background = 1, 2
+markers[coins < 30] = background
+markers[coins > 150] = foreground
+
+ws = watershed(edges, markers)
+seg1 = nd.label(ws == foreground)[0]
+
+# make segmentation using SLIC superpixels
+
+# make the RGB equivalent of `coins`
+coins_colour = np.tile(coins[..., np.newaxis], (1, 1, 3))
+seg2 = slic(coins_colour, n_segments=30, max_iter=160, sigma=1, ratio=9,
+            convert2lab=False)
+
+# combine the two
+segj = join_segmentations(seg1, seg2)
+
+### Display the result ###
+
+# make a random colormap for a set number of values
+def random_cmap(im):
+    np.random.seed(9)
+    cmap_array = np.concatenate(
+        (np.zeros((1, 3)), np.random.rand(np.ceil(im.max()), 3)))
+    return mpl.colors.ListedColormap(cmap_array)
+
+# show the segmentations
+fig, axes = plt.subplots(ncols=4, figsize=(9, 2.5))
+axes[0].imshow(coins, cmap=plt.cm.gray, interpolation='nearest')
+axes[0].set_title('Image')
+axes[1].imshow(seg1, cmap=random_cmap(seg1), interpolation='nearest')
+axes[1].set_title('Sobel+Watershed')
+axes[2].imshow(seg2, cmap=random_cmap(seg2), interpolation='nearest')
+axes[2].set_title('SLIC superpixels')
+axes[3].imshow(segj, cmap=random_cmap(segj), interpolation='nearest')
+axes[3].set_title('Join')
+
+for ax in axes:
+    ax.axis('off')
+plt.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0, right=1)
+plt.show()
diff --git a/skimage/segmentation/__init__.py b/skimage/segmentation/__init__.py
index 7bd37026..a0fab77a 100644
--- a/skimage/segmentation/__init__.py
+++ b/skimage/segmentation/__init__.py
@@ -4,3 +4,4 @@ from ._slic import slic
 from ._quickshift import quickshift
 from .boundaries import find_boundaries, visualize_boundaries, mark_boundaries
 from ._clear_border import clear_border
+from ._join import join_segmentations, relabel_from_one
diff --git a/skimage/segmentation/_join.py b/skimage/segmentation/_join.py
new file mode 100644
index 00000000..454da71e
--- /dev/null
+++ b/skimage/segmentation/_join.py
@@ -0,0 +1,93 @@
+import numpy as np
+
+def join_segmentations(s1, s2):
+    """Return the join of the two input segmentations.
+
+    The join J of S1 and S2 is defined as the segmentation in which two voxels
+    are in the same segment in J if and only if they are in the same segment
+    in *both* S1 and S2.
+
+    Parameters
+    ----------
+    s1, s2 : numpy arrays
+        s1 and s2 are label fields of the same shape.
+
+    Returns
+    -------
+    j : numpy array
+        The join segmentation of s1 and s2.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> from skimage.segmentation import join_segmentations
+    >>> s1 = np.array([[0, 0, 1, 1],
+    ...                [0, 2, 1, 1],
+    ...                [2, 2, 2, 1]])
+    >>> s2 = np.array([[0, 1, 1, 0],
+    ...                [0, 1, 1, 0],
+    ...                [0, 1, 1, 1]])
+    >>> join_segmentations(s1, s2)
+    array([[0, 1, 3, 2],
+           [0, 5, 3, 2],
+           [4, 5, 5, 3]])
+    """
+    if s1.shape != s2.shape:
+        raise ValueError("Cannot join segmentations of different shape. " +
+                         "s1.shape: %s, s2.shape: %s" % (s1.shape, s2.shape))
+    s1 = relabel_from_one(s1)[0]
+    s2 = relabel_from_one(s2)[0]
+    j = (s2.max() + 1) * s1 + s2
+    j = relabel_from_one(j)[0]
+    return j
+
+def relabel_from_one(label_field):
+    """Convert labels in an arbitrary label field to {1, ... number_of_labels}.
+
+    This function also returns the forward map (mapping the original labels to
+    the reduced labels) and the inverse map (mapping the reduced labels back
+    to the original ones).
+
+    Parameters
+    ----------
+    label_field : numpy ndarray (integer type)
+
+    Returns
+    -------
+    relabeled : numpy array of same shape as ar
+    forward_map : 1d numpy array of length np.unique(ar) + 1
+    inverse_map : 1d numpy array of length len(np.unique(ar))
+        The length is len(np.unique(ar)) + 1 if 0 is not in np.unique(ar)
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> from skimage.segmentation import relabel_from_one
+    >>> label_field = array([1, 1, 5, 5, 8, 99, 42])
+    >>> relab, fw, inv = relabel_from_one(label_field)
+    >>> relab
+    array([1, 1, 2, 2, 3, 5, 4])
+    >>> fw
+    array([0, 1, 0, 0, 0, 2, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 0, 0,
+           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+           0, 0, 0, 0, 0, 0, 0, 5])
+    >>> inv
+    array([ 0,  1,  5,  8, 42, 99])
+    >>> (fw[label_field] == relab).all()
+    True
+    >>> (inv[relab] == label_field).all()
+    True
+    """
+    labels = np.unique(label_field)
+    labels0 = labels[labels != 0]
+    m = labels.max()
+    if m == len(labels0): # nothing to do, already 1...n labels
+        return label_field, labels, labels
+    forward_map = np.zeros(m+1, int)
+    forward_map[labels0] = np.arange(1, len(labels0) + 1)
+    if not (labels == 0).any():
+        labels = np.concatenate(([0], labels))
+    inverse_map = labels
+    return forward_map[label_field], forward_map, inverse_map
diff --git a/skimage/segmentation/tests/test_join.py b/skimage/segmentation/tests/test_join.py
new file mode 100644
index 00000000..f03244e9
--- /dev/null
+++ b/skimage/segmentation/tests/test_join.py
@@ -0,0 +1,40 @@
+import numpy as np
+from numpy.testing import assert_array_equal, assert_raises
+from skimage.segmentation import join_segmentations, relabel_from_one
+
+def test_join_segmentations():
+    s1 = np.array([[0, 0, 1, 1],
+                   [0, 2, 1, 1],
+                   [2, 2, 2, 1]])
+    s2 = np.array([[0, 1, 1, 0],
+                   [0, 1, 1, 0],
+                   [0, 1, 1, 1]])
+
+    # test correct join
+    # NOTE: technically, equality to j_ref is not required, only that there
+    # is a one-to-one mapping between j and j_ref. I don't know of an easy way
+    # to check this (i.e. not as error-prone as the function being tested)
+    j = join_segmentations(s1, s2)
+    j_ref = np.array([[0, 1, 3, 2],
+                      [0, 5, 3, 2],
+                      [4, 5, 5, 3]])
+    assert_array_equal(j, j_ref)
+
+    # test correct exception when arrays are different shapes
+    s3 = np.array([[0, 0, 1, 1], [0, 2, 2, 1]])
+    assert_raises(ValueError, join_segmentations, s1, s3)
+
+def test_relabel_from_one():
+    ar = np.array([1, 1, 5, 5, 8, 99, 42])
+    ar_relab, fw, inv = relabel_from_one(ar)
+    ar_relab_ref = np.array([1, 1, 2, 2, 3, 5, 4])
+    assert_array_equal(ar_relab, ar_relab_ref)
+    fw_ref = np.zeros(100, int)
+    fw_ref[1] = 1; fw_ref[5] = 2; fw_ref[8] = 3; fw_ref[42] = 4; fw_ref[99] = 5
+    assert_array_equal(fw, fw_ref)
+    inv_ref = np.array([0,  1,  5,  8, 42, 99])
+    assert_array_equal(inv, inv_ref)
+
+
+if __name__ == "__main__":
+    np.testing.run_module_suite()