_join.py 7.0 KB

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