1import numpy as np
2from ..util._map_array import map_array, ArrayMap
3
4
5def join_segmentations(s1, s2):
6    """Return the join of the two input segmentations.
7
8    The join J of S1 and S2 is defined as the segmentation in which two
9    voxels are in the same segment if and only if they are in the same
10    segment in *both* S1 and S2.
11
12    Parameters
13    ----------
14    s1, s2 : numpy arrays
15        s1 and s2 are label fields of the same shape.
16
17    Returns
18    -------
19    j : numpy array
20        The join segmentation of s1 and s2.
21
22    Examples
23    --------
24    >>> from skimage.segmentation import join_segmentations
25    >>> s1 = np.array([[0, 0, 1, 1],
26    ...                [0, 2, 1, 1],
27    ...                [2, 2, 2, 1]])
28    >>> s2 = np.array([[0, 1, 1, 0],
29    ...                [0, 1, 1, 0],
30    ...                [0, 1, 1, 1]])
31    >>> join_segmentations(s1, s2)
32    array([[0, 1, 3, 2],
33           [0, 5, 3, 2],
34           [4, 5, 5, 3]])
35    """
36    if s1.shape != s2.shape:
37        raise ValueError("Cannot join segmentations of different shape. " +
38                         "s1.shape: %s, s2.shape: %s" % (s1.shape, s2.shape))
39    s1 = relabel_sequential(s1)[0]
40    s2 = relabel_sequential(s2)[0]
41    j = (s2.max() + 1) * s1 + s2
42    j = relabel_sequential(j)[0]
43    return j
44
45
46def relabel_sequential(label_field, offset=1):
47    """Relabel arbitrary labels to {`offset`, ... `offset` + number_of_labels}.
48
49    This function also returns the forward map (mapping the original labels to
50    the reduced labels) and the inverse map (mapping the reduced labels back
51    to the original ones).
52
53    Parameters
54    ----------
55    label_field : numpy array of int, arbitrary shape
56        An array of labels, which must be non-negative integers.
57    offset : int, optional
58        The return labels will start at `offset`, which should be
59        strictly positive.
60
61    Returns
62    -------
63    relabeled : numpy array of int, same shape as `label_field`
64        The input label field with labels mapped to
65        {offset, ..., number_of_labels + offset - 1}.
66        The data type will be the same as `label_field`, except when
67        offset + number_of_labels causes overflow of the current data type.
68    forward_map : ArrayMap
69        The map from the original label space to the returned label
70        space. Can be used to re-apply the same mapping. See examples
71        for usage. The output data type will be the same as `relabeled`.
72    inverse_map : ArrayMap
73        The map from the new label space to the original space. This
74        can be used to reconstruct the original label field from the
75        relabeled one. The output data type will be the same as `label_field`.
76
77    Notes
78    -----
79    The label 0 is assumed to denote the background and is never remapped.
80
81    The forward map can be extremely big for some inputs, since its
82    length is given by the maximum of the label field. However, in most
83    situations, ``label_field.max()`` is much smaller than
84    ``label_field.size``, and in these cases the forward map is
85    guaranteed to be smaller than either the input or output images.
86
87    Examples
88    --------
89    >>> from skimage.segmentation import relabel_sequential
90    >>> label_field = np.array([1, 1, 5, 5, 8, 99, 42])
91    >>> relab, fw, inv = relabel_sequential(label_field)
92    >>> relab
93    array([1, 1, 2, 2, 3, 5, 4])
94    >>> print(fw)
95    ArrayMap:
96      1 → 1
97      5 → 2
98      8 → 3
99      42 → 4
100      99 → 5
101    >>> np.array(fw)
102    array([0, 1, 0, 0, 0, 2, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
103           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0,
104           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
105           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
106           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5])
107    >>> np.array(inv)
108    array([ 0,  1,  5,  8, 42, 99])
109    >>> (fw[label_field] == relab).all()
110    True
111    >>> (inv[relab] == label_field).all()
112    True
113    >>> relab, fw, inv = relabel_sequential(label_field, offset=5)
114    >>> relab
115    array([5, 5, 6, 6, 7, 9, 8])
116    """
117    if offset <= 0:
118        raise ValueError("Offset must be strictly positive.")
119    if np.min(label_field) < 0:
120        raise ValueError("Cannot relabel array that contains negative values.")
121    offset = int(offset)
122    in_vals = np.unique(label_field)
123    if in_vals[0] == 0:
124        # always map 0 to 0
125        out_vals = np.concatenate(
126            [[0], np.arange(offset, offset+len(in_vals)-1)]
127        )
128    else:
129        out_vals = np.arange(offset, offset+len(in_vals))
130    input_type = label_field.dtype
131
132    # Some logic to determine the output type:
133    #  - we don't want to return a smaller output type than the input type,
134    #  ie if we get uint32 as labels input, don't return a uint8 array.
135    #  - but, in some cases, using the input type could result in overflow. The
136    #  input type could be a signed integer (e.g. int32) but
137    #  `np.min_scalar_type` will always return an unsigned type. We check for
138    #  that by casting the largest output value to the input type. If it is
139    #  unchanged, we use the input type, else we use the unsigned minimum
140    #  required type
141    required_type = np.min_scalar_type(out_vals[-1])
142    if input_type.itemsize < required_type.itemsize:
143        output_type = required_type
144    else:
145        if input_type.type(out_vals[-1]) == out_vals[-1]:
146            output_type = input_type
147        else:
148            output_type = required_type
149    out_array = np.empty(label_field.shape, dtype=output_type)
150    out_vals = out_vals.astype(output_type)
151    map_array(label_field, in_vals, out_vals, out=out_array)
152    fw_map = ArrayMap(in_vals, out_vals)
153    inv_map = ArrayMap(out_vals, in_vals)
154    return out_array, fw_map, inv_map
155