From dc377970f520f8f04a2b3d879983c3a1171518e7 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <jni@janelia.hhmi.org>
Date: Wed, 28 Nov 2012 12:22:14 +1100
Subject: [PATCH] Add functions to join and relabel segmentations

The join of two segmentations is the segmentation in which two voxels
are in the same segment if and only if they are in the same segment in
both input segmentations.
---
 skimage/segmentation/__init__.py |  1 +
 skimage/segmentation/_join.py    | 93 ++++++++++++++++++++++++++++++++
 2 files changed, 94 insertions(+)
 create mode 100644 skimage/segmentation/_join.py

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..5cd25b80
--- /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(ar):
+    """Convert array ar of arbitrary labels to labels 1...len(np.unique(ar))+1
+
+    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
+    ----------
+    ar : numpy ndarray (integer type)
+
+    Returns
+    -------
+    ar_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
+    >>> ar = array([1, 1, 5, 5, 8, 99, 42])
+    >>> ar_relab, fw, inv = relabel_from_one(ar)
+    >>> ar_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[ar] == ar_relab).all()
+    True
+    >>> (inv[ar_relab] == ar).all()
+    True
+    """
+    labels = np.unique(ar)
+    labels0 = labels[labels != 0]
+    m = labels.max()
+    if m == len(labels0): # nothing to do, already 1...n labels
+        return ar, 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[ar], forward_map, inverse_map