boundaries.py 9.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240
  1. import numpy as np
  2. from scipy import ndimage as ndi
  3. from .._shared.utils import _supported_float_type
  4. from ..morphology import dilation, erosion, footprint_rectangle
  5. from ..util import img_as_float, view_as_windows
  6. from ..color import gray2rgb
  7. def _find_boundaries_subpixel(label_img):
  8. """See ``find_boundaries(..., mode='subpixel')``.
  9. Notes
  10. -----
  11. This function puts in an empty row and column between each *actual*
  12. row and column of the image, for a corresponding shape of ``2s - 1``
  13. for every image dimension of size ``s``. These "interstitial" rows
  14. and columns are filled as ``True`` if they separate two labels in
  15. `label_img`, ``False`` otherwise.
  16. I used ``view_as_windows`` to get the neighborhood of each pixel.
  17. Then I check whether there are two labels or more in that
  18. neighborhood.
  19. """
  20. ndim = label_img.ndim
  21. max_label = np.iinfo(label_img.dtype).max
  22. label_img_expanded = np.zeros(
  23. [(2 * s - 1) for s in label_img.shape], label_img.dtype
  24. )
  25. pixels = (slice(None, None, 2),) * ndim
  26. label_img_expanded[pixels] = label_img
  27. edges = np.ones(label_img_expanded.shape, dtype=bool)
  28. edges[pixels] = False
  29. label_img_expanded[edges] = max_label
  30. windows = view_as_windows(np.pad(label_img_expanded, 1, mode='edge'), (3,) * ndim)
  31. boundaries = np.zeros_like(edges)
  32. for index in np.ndindex(label_img_expanded.shape):
  33. if edges[index]:
  34. values = np.unique(windows[index].ravel())
  35. if len(values) > 2: # single value and max_label
  36. boundaries[index] = True
  37. return boundaries
  38. def find_boundaries(label_img, connectivity=1, mode='thick', background=0):
  39. """Return bool array where boundaries between labeled regions are True.
  40. Parameters
  41. ----------
  42. label_img : array of int or bool
  43. An array in which different regions are labeled with either different
  44. integers or boolean values.
  45. connectivity : int in {1, ..., `label_img.ndim`}, optional
  46. A pixel is considered a boundary pixel if any of its neighbors
  47. has a different label. `connectivity` controls which pixels are
  48. considered neighbors. A connectivity of 1 (default) means
  49. pixels sharing an edge (in 2D) or a face (in 3D) will be
  50. considered neighbors. A connectivity of `label_img.ndim` means
  51. pixels sharing a corner will be considered neighbors.
  52. mode : string in {'thick', 'inner', 'outer', 'subpixel'}
  53. How to mark the boundaries:
  54. - thick: any pixel not completely surrounded by pixels of the
  55. same label (defined by `connectivity`) is marked as a boundary.
  56. This results in boundaries that are 2 pixels thick.
  57. - inner: outline the pixels *just inside* of objects, leaving
  58. background pixels untouched.
  59. - outer: outline pixels in the background around object
  60. boundaries. When two objects touch, their boundary is also
  61. marked.
  62. - subpixel: return a doubled image, with pixels *between* the
  63. original pixels marked as boundary where appropriate.
  64. background : int, optional
  65. For modes 'inner' and 'outer', a definition of a background
  66. label is required. See `mode` for descriptions of these two.
  67. Returns
  68. -------
  69. boundaries : array of bool, same shape as `label_img`
  70. A bool image where ``True`` represents a boundary pixel. For
  71. `mode` equal to 'subpixel', ``boundaries.shape[i]`` is equal
  72. to ``2 * label_img.shape[i] - 1`` for all ``i`` (a pixel is
  73. inserted in between all other pairs of pixels).
  74. Examples
  75. --------
  76. >>> labels = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  77. ... [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  78. ... [0, 0, 0, 0, 0, 5, 5, 5, 0, 0],
  79. ... [0, 0, 1, 1, 1, 5, 5, 5, 0, 0],
  80. ... [0, 0, 1, 1, 1, 5, 5, 5, 0, 0],
  81. ... [0, 0, 1, 1, 1, 5, 5, 5, 0, 0],
  82. ... [0, 0, 0, 0, 0, 5, 5, 5, 0, 0],
  83. ... [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  84. ... [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=np.uint8)
  85. >>> find_boundaries(labels, mode='thick').astype(np.uint8)
  86. array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  87. [0, 0, 0, 0, 0, 1, 1, 1, 0, 0],
  88. [0, 0, 1, 1, 1, 1, 1, 1, 1, 0],
  89. [0, 1, 1, 1, 1, 1, 0, 1, 1, 0],
  90. [0, 1, 1, 0, 1, 1, 0, 1, 1, 0],
  91. [0, 1, 1, 1, 1, 1, 0, 1, 1, 0],
  92. [0, 0, 1, 1, 1, 1, 1, 1, 1, 0],
  93. [0, 0, 0, 0, 0, 1, 1, 1, 0, 0],
  94. [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=uint8)
  95. >>> find_boundaries(labels, mode='inner').astype(np.uint8)
  96. array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  97. [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  98. [0, 0, 0, 0, 0, 1, 1, 1, 0, 0],
  99. [0, 0, 1, 1, 1, 1, 0, 1, 0, 0],
  100. [0, 0, 1, 0, 1, 1, 0, 1, 0, 0],
  101. [0, 0, 1, 1, 1, 1, 0, 1, 0, 0],
  102. [0, 0, 0, 0, 0, 1, 1, 1, 0, 0],
  103. [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  104. [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=uint8)
  105. >>> find_boundaries(labels, mode='outer').astype(np.uint8)
  106. array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  107. [0, 0, 0, 0, 0, 1, 1, 1, 0, 0],
  108. [0, 0, 1, 1, 1, 1, 0, 0, 1, 0],
  109. [0, 1, 0, 0, 1, 1, 0, 0, 1, 0],
  110. [0, 1, 0, 0, 1, 1, 0, 0, 1, 0],
  111. [0, 1, 0, 0, 1, 1, 0, 0, 1, 0],
  112. [0, 0, 1, 1, 1, 1, 0, 0, 1, 0],
  113. [0, 0, 0, 0, 0, 1, 1, 1, 0, 0],
  114. [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=uint8)
  115. >>> labels_small = labels[::2, ::3]
  116. >>> labels_small
  117. array([[0, 0, 0, 0],
  118. [0, 0, 5, 0],
  119. [0, 1, 5, 0],
  120. [0, 0, 5, 0],
  121. [0, 0, 0, 0]], dtype=uint8)
  122. >>> find_boundaries(labels_small, mode='subpixel').astype(np.uint8)
  123. array([[0, 0, 0, 0, 0, 0, 0],
  124. [0, 0, 0, 1, 1, 1, 0],
  125. [0, 0, 0, 1, 0, 1, 0],
  126. [0, 1, 1, 1, 0, 1, 0],
  127. [0, 1, 0, 1, 0, 1, 0],
  128. [0, 1, 1, 1, 0, 1, 0],
  129. [0, 0, 0, 1, 0, 1, 0],
  130. [0, 0, 0, 1, 1, 1, 0],
  131. [0, 0, 0, 0, 0, 0, 0]], dtype=uint8)
  132. >>> bool_image = np.array([[False, False, False, False, False],
  133. ... [False, False, False, False, False],
  134. ... [False, False, True, True, True],
  135. ... [False, False, True, True, True],
  136. ... [False, False, True, True, True]],
  137. ... dtype=bool)
  138. >>> find_boundaries(bool_image)
  139. array([[False, False, False, False, False],
  140. [False, False, True, True, True],
  141. [False, True, True, True, True],
  142. [False, True, True, False, False],
  143. [False, True, True, False, False]])
  144. """
  145. if label_img.dtype == 'bool':
  146. label_img = label_img.astype(np.uint8)
  147. ndim = label_img.ndim
  148. footprint = ndi.generate_binary_structure(ndim, connectivity)
  149. if mode != 'subpixel':
  150. boundaries = dilation(label_img, footprint) != erosion(label_img, footprint)
  151. if mode == 'inner':
  152. foreground_image = label_img != background
  153. boundaries &= foreground_image
  154. elif mode == 'outer':
  155. max_label = np.iinfo(label_img.dtype).max
  156. background_image = label_img == background
  157. footprint = ndi.generate_binary_structure(ndim, ndim)
  158. inverted_background = np.array(label_img, copy=True)
  159. inverted_background[background_image] = max_label
  160. adjacent_objects = (
  161. dilation(label_img, footprint)
  162. != erosion(inverted_background, footprint)
  163. ) & ~background_image
  164. boundaries &= background_image | adjacent_objects
  165. return boundaries
  166. else:
  167. boundaries = _find_boundaries_subpixel(label_img)
  168. return boundaries
  169. def mark_boundaries(
  170. image,
  171. label_img,
  172. color=(1, 1, 0),
  173. outline_color=None,
  174. mode='outer',
  175. background_label=0,
  176. ):
  177. """Return image with boundaries between labeled regions highlighted.
  178. Parameters
  179. ----------
  180. image : (M, N[, 3]) array
  181. Grayscale or RGB image.
  182. label_img : (M, N) array of int
  183. Label array where regions are marked by different integer values.
  184. color : length-3 sequence, optional
  185. RGB color of boundaries in the output image.
  186. outline_color : length-3 sequence, optional
  187. RGB color surrounding boundaries in the output image. If None, no
  188. outline is drawn.
  189. mode : string in {'thick', 'inner', 'outer', 'subpixel'}, optional
  190. The mode for finding boundaries.
  191. background_label : int, optional
  192. Which label to consider background (this is only useful for
  193. modes ``inner`` and ``outer``).
  194. Returns
  195. -------
  196. marked : (M, N, 3) array of float
  197. An image in which the boundaries between labels are
  198. superimposed on the original image.
  199. See Also
  200. --------
  201. find_boundaries
  202. """
  203. float_dtype = _supported_float_type(image.dtype)
  204. marked = img_as_float(image, force_copy=True)
  205. marked = marked.astype(float_dtype, copy=False)
  206. if marked.ndim == 2:
  207. marked = gray2rgb(marked)
  208. if mode == 'subpixel':
  209. # Here, we want to interpose an extra line of pixels between
  210. # each original line - except for the last axis which holds
  211. # the RGB information. ``ndi.zoom`` then performs the (cubic)
  212. # interpolation, filling in the values of the interposed pixels
  213. marked = ndi.zoom(
  214. marked, [2 - 1 / s for s in marked.shape[:-1]] + [1], mode='mirror'
  215. )
  216. boundaries = find_boundaries(label_img, mode=mode, background=background_label)
  217. if outline_color is not None:
  218. outlines = dilation(boundaries, footprint_rectangle((3, 3)))
  219. marked[outlines] = outline_color
  220. marked[boundaries] = color
  221. return marked