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