misc.py 18 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509
  1. """Miscellaneous morphology functions."""
  2. import numpy as np
  3. import functools
  4. import warnings
  5. from scipy import ndimage as ndi
  6. from scipy.spatial import cKDTree
  7. from .._shared.utils import warn, deprecate_parameter, DEPRECATED
  8. from ._misc_cy import _remove_objects_by_distance
  9. # Our function names don't exactly correspond to ndimages.
  10. # This dictionary translates from our names to scipy's.
  11. funcs = ('erosion', 'dilation', 'opening', 'closing')
  12. skimage2ndimage = {x: 'grey_' + x for x in funcs}
  13. # These function names are the same in ndimage.
  14. funcs = (
  15. 'black_tophat',
  16. 'white_tophat',
  17. )
  18. skimage2ndimage.update({x: x for x in funcs})
  19. def default_footprint(func):
  20. """Decorator to add a default footprint to morphology functions.
  21. Parameters
  22. ----------
  23. func : function
  24. A morphology function such as erosion, dilation, opening, closing,
  25. white_tophat, or black_tophat.
  26. Returns
  27. -------
  28. func_out : function
  29. The function, using a default footprint of same dimension
  30. as the input image with connectivity 1.
  31. """
  32. @functools.wraps(func)
  33. def func_out(image, footprint=None, *args, **kwargs):
  34. if footprint is None:
  35. footprint = ndi.generate_binary_structure(image.ndim, 1)
  36. return func(image, footprint=footprint, *args, **kwargs)
  37. return func_out
  38. def _check_dtype_supported(ar):
  39. # Should use `issubdtype` for bool below, but there's a bug in numpy 1.7
  40. if not (ar.dtype == bool or np.issubdtype(ar.dtype, np.integer)):
  41. raise TypeError(
  42. "Only bool or integer image types are supported. " f"Got {ar.dtype}."
  43. )
  44. @deprecate_parameter(
  45. deprecated_name="min_size",
  46. new_name="max_size",
  47. start_version="0.26.0",
  48. stop_version="2.0.0",
  49. template=f"{deprecate_parameter.replace_parameter_template} "
  50. "Note that the new threshold removes objects smaller than **or equal to** "
  51. "its value, while the previous parameter only removed smaller ones.",
  52. )
  53. def remove_small_objects(
  54. ar, min_size=DEPRECATED, connectivity=1, *, max_size=64, out=None
  55. ):
  56. """Remove objects smaller than the specified size.
  57. Expects `ar` to be an array with labeled objects, and removes objects
  58. smaller than or equal to `max_size`. If `ar` is bool, the image is first
  59. labeled. This leads to potentially different behavior for bool vs. 0-and-1
  60. arrays.
  61. Parameters
  62. ----------
  63. ar : ndarray (arbitrary shape, int or bool type)
  64. The array containing the objects of interest. If the array type is
  65. int, the ints must be non-negative.
  66. max_size : int, optional (default: 64)
  67. Remove objects whose contiguous area (or volume, in N-D) contains this
  68. number of pixels or fewer.
  69. .. versionadded:: 0.26
  70. To make the naming clearer, replaces deprecated `min_size`
  71. which only removed objects strictly smaller than its size.
  72. connectivity : int, {1, 2, ..., ar.ndim}, optional (default: 1)
  73. The connectivity defining the neighborhood of a pixel. Used during
  74. labelling if `ar` is bool.
  75. out : ndarray
  76. Array of the same shape as `ar`, into which the output is
  77. placed. By default, a new array is created.
  78. Raises
  79. ------
  80. TypeError
  81. If the input array is of an invalid type, such as float or string.
  82. ValueError
  83. If the input array contains negative values.
  84. Returns
  85. -------
  86. out : ndarray, same shape and type as input `ar`
  87. The input array with small connected components removed.
  88. See Also
  89. --------
  90. skimage.morphology.remove_small_holes
  91. skimage.morphology.remove_objects_by_distance
  92. Examples
  93. --------
  94. >>> from skimage import morphology
  95. >>> a = np.array([[0, 0, 0, 1, 0],
  96. ... [1, 1, 1, 0, 0],
  97. ... [1, 1, 1, 0, 1]], bool)
  98. >>> b = morphology.remove_small_objects(a, max_size=5)
  99. >>> b
  100. array([[False, False, False, False, False],
  101. [ True, True, True, False, False],
  102. [ True, True, True, False, False]])
  103. >>> c = morphology.remove_small_objects(a, max_size=6, connectivity=2)
  104. >>> c
  105. array([[False, False, False, True, False],
  106. [ True, True, True, False, False],
  107. [ True, True, True, False, False]])
  108. >>> d = morphology.remove_small_objects(a, max_size=5, out=a)
  109. >>> d is a
  110. True
  111. """
  112. # Raising type error if not int or bool
  113. _check_dtype_supported(ar)
  114. if out is None:
  115. out = ar.copy()
  116. else:
  117. out[:] = ar
  118. if max_size == 0: # shortcut for efficiency
  119. return out
  120. if out.dtype == bool:
  121. footprint = ndi.generate_binary_structure(ar.ndim, connectivity)
  122. ccs = np.zeros_like(ar, dtype=np.int32)
  123. ndi.label(ar, footprint, output=ccs)
  124. else:
  125. ccs = out
  126. try:
  127. component_sizes = np.bincount(ccs.ravel())
  128. except ValueError:
  129. raise ValueError(
  130. "Negative value labels are not supported. Try "
  131. "relabeling the input with `scipy.ndimage.label` or "
  132. "`skimage.morphology.label`."
  133. )
  134. if len(component_sizes) == 2 and out.dtype != bool:
  135. warn(
  136. "Only one label was provided to `remove_small_objects`. "
  137. "Did you mean to use a boolean array?"
  138. )
  139. if min_size is not DEPRECATED:
  140. # Exclusive threshold is deprecated behavior
  141. too_small = component_sizes < min_size
  142. else:
  143. # New behavior uses inclusive threshold
  144. too_small = component_sizes <= max_size
  145. too_small_mask = too_small[ccs]
  146. out[too_small_mask] = 0
  147. return out
  148. @deprecate_parameter(
  149. deprecated_name="area_threshold",
  150. new_name="max_size",
  151. start_version="0.26.0",
  152. stop_version="2.0.0",
  153. template=f"{deprecate_parameter.replace_parameter_template} "
  154. "Note that the new threshold removes objects smaller than **or equal to** "
  155. "its value, while the previous parameter only removed smaller ones.",
  156. )
  157. def remove_small_holes(
  158. ar, area_threshold=DEPRECATED, connectivity=1, *, max_size=64, out=None
  159. ):
  160. """Remove contiguous holes smaller than the specified size.
  161. Parameters
  162. ----------
  163. ar : ndarray (arbitrary shape, int or bool type)
  164. The array containing the connected components of interest.
  165. max_size : int, optional (default: 64)
  166. Remove holes whose contiguous area (or volume, in N-D) contains this
  167. number of pixels or fewer.
  168. .. versionadded:: 0.26
  169. To make the naming clearer, replaces deprecated `area_threshold`
  170. which only removed holes strictly smaller than its size.
  171. connectivity : int, {1, 2, ..., ar.ndim}, optional (default: 1)
  172. The connectivity defining the neighborhood of a pixel.
  173. out : ndarray
  174. Array of the same shape as `ar` and bool dtype, into which the
  175. output is placed. By default, a new array is created.
  176. Raises
  177. ------
  178. TypeError
  179. If the input array is of an invalid type, such as float or string.
  180. ValueError
  181. If the input array contains negative values.
  182. Returns
  183. -------
  184. out : ndarray, same shape and type as input `ar`
  185. The input array with small holes within connected components removed.
  186. See Also
  187. --------
  188. skimage.morphology.remove_small_objects
  189. skimage.morphology.remove_objects_by_distance
  190. Examples
  191. --------
  192. >>> from skimage import morphology
  193. >>> a = np.array([[1, 1, 1, 1, 1, 0],
  194. ... [1, 1, 1, 0, 1, 0],
  195. ... [1, 0, 0, 1, 1, 0],
  196. ... [1, 1, 1, 1, 1, 0]], bool)
  197. >>> b = morphology.remove_small_holes(a, max_size=1)
  198. >>> b
  199. array([[ True, True, True, True, True, False],
  200. [ True, True, True, True, True, False],
  201. [ True, False, False, True, True, False],
  202. [ True, True, True, True, True, False]])
  203. >>> c = morphology.remove_small_holes(a, max_size=1, connectivity=2)
  204. >>> c
  205. array([[ True, True, True, True, True, False],
  206. [ True, True, True, False, True, False],
  207. [ True, False, False, True, True, False],
  208. [ True, True, True, True, True, False]])
  209. >>> d = morphology.remove_small_holes(a, max_size=1, out=a)
  210. >>> d is a
  211. True
  212. Notes
  213. -----
  214. If the array type is int, it is assumed that it contains already-labeled
  215. objects. The labels are not kept in the output image (this function always
  216. outputs a bool image). It is suggested that labeling is completed after
  217. using this function.
  218. """
  219. _check_dtype_supported(ar)
  220. # Creates warning if image is an integer image
  221. if ar.dtype != bool:
  222. warn(
  223. "Any labeled images will be returned as a boolean array. "
  224. "Did you mean to use a boolean array?",
  225. UserWarning,
  226. )
  227. if out is not None:
  228. if out.dtype != bool:
  229. raise TypeError("out dtype must be bool")
  230. else:
  231. out = ar.astype(bool, copy=True)
  232. # Creating the inverse of ar
  233. np.logical_not(ar, out=out)
  234. # removing small objects from the inverse of ar
  235. with warnings.catch_warnings():
  236. warnings.filterwarnings(
  237. "ignore",
  238. message="Parameter `min_size` is deprecated",
  239. category=FutureWarning,
  240. )
  241. out = remove_small_objects(
  242. out,
  243. min_size=area_threshold,
  244. max_size=max_size,
  245. connectivity=connectivity,
  246. out=out,
  247. )
  248. np.logical_not(out, out=out)
  249. return out
  250. def remove_objects_by_distance(
  251. label_image,
  252. min_distance,
  253. *,
  254. priority=None,
  255. p_norm=2,
  256. spacing=None,
  257. out=None,
  258. ):
  259. """Remove objects, in specified order, until remaining are a minimum distance apart.
  260. Remove labeled objects from an image until the remaining ones are spaced
  261. more than a given distance from one another. By default, smaller objects
  262. are removed first.
  263. Parameters
  264. ----------
  265. label_image : ndarray of integers
  266. An n-dimensional array containing object labels, e.g. as returned by
  267. :func:`~.label`. A value of zero is considered background, all other
  268. object IDs must be positive integers.
  269. min_distance : int or float
  270. Remove objects whose distance to other objects is not greater than this
  271. positive value. Objects with a lower `priority` are removed first.
  272. priority : ndarray, optional
  273. Defines the priority with which objects are removed. Expects a
  274. 1-dimensional array of length
  275. :func:`np.amax(label_image) + 1 <numpy.amax>` that contains the priority
  276. for each object's label at the respective index. Objects with a lower value
  277. are removed first until all remaining objects fulfill the distance
  278. requirement. If not given, priority is given to objects with a higher
  279. number of samples and their label value second.
  280. p_norm : int or float, optional
  281. The Minkowski distance of order p, used to calculate the distance
  282. between objects. The default ``2`` corresponds to the Euclidean
  283. distance, ``1`` to the "Manhattan" distance, and ``np.inf`` to the
  284. Chebyshev distance.
  285. spacing : sequence of float, optional
  286. The pixel spacing along each axis of `label_image`. If not specified,
  287. a grid spacing of unity (1) is implied.
  288. out : ndarray, optional
  289. Array of the same shape and dtype as `image`, into which the output is
  290. placed. By default, a new array is created.
  291. Returns
  292. -------
  293. out : ndarray
  294. Array of the same shape as `label_image`, for which objects that violate
  295. the `min_distance` condition were removed.
  296. See Also
  297. --------
  298. skimage.morphology.remove_small_objects
  299. Remove objects smaller than the specified size.
  300. skimage.morphology.remove_small_holes
  301. Remove holes smaller than the specified size.
  302. Notes
  303. -----
  304. The basic steps of this algorithm work as follows:
  305. 1. Find the indices for of all given objects and separate them depending on
  306. if they point to an object's border or not.
  307. 2. Sort indices by their label value, ensuring that indices which point to
  308. the same object are next to each other. This optimization allows finding
  309. all parts of an object, simply by stepping to the neighboring indices.
  310. 3. Sort boundary indices by `priority`. Use a stable-sort to preserve the
  311. ordering from the previous sorting step. If `priority` is not given,
  312. use :func:`numpy.bincount` as a fallback.
  313. 4. Construct a :class:`scipy.spatial.cKDTree` from the boundary indices.
  314. 5. Iterate across boundary indices in priority-sorted order, and query the
  315. kd-tree for objects that are too close. Remove ones that are and don't
  316. take them into account when evaluating other objects later on.
  317. The performance of this algorithm depends on the number of samples in
  318. `label_image` that belong to an object's border.
  319. Examples
  320. --------
  321. >>> import skimage as ski
  322. >>> ski.morphology.remove_objects_by_distance(np.array([2, 0, 1, 1]), 2)
  323. array([0, 0, 1, 1])
  324. >>> ski.morphology.remove_objects_by_distance(
  325. ... np.array([2, 0, 1, 1]), 2, priority=np.array([0, 1, 9])
  326. ... )
  327. array([2, 0, 0, 0])
  328. >>> label_image = np.array(
  329. ... [[8, 0, 0, 0, 0, 0, 0, 0, 0, 9, 9],
  330. ... [8, 8, 8, 0, 0, 0, 0, 0, 0, 9, 9],
  331. ... [0, 0, 0, 0, 0, 0, 0, 0, 9, 0, 0],
  332. ... [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  333. ... [0, 0, 3, 0, 0, 0, 1, 0, 0, 0, 0],
  334. ... [2, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
  335. ... [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
  336. ... [0, 0, 0, 0, 0, 0, 0, 0, 0, 7, 7]]
  337. ... )
  338. >>> ski.morphology.remove_objects_by_distance(
  339. ... label_image, min_distance=3
  340. ... )
  341. array([[8, 0, 0, 0, 0, 0, 0, 0, 0, 9, 9],
  342. [8, 8, 8, 0, 0, 0, 0, 0, 0, 9, 9],
  343. [0, 0, 0, 0, 0, 0, 0, 0, 9, 0, 0],
  344. [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  345. [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  346. [2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  347. [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  348. [0, 0, 0, 0, 0, 0, 0, 0, 0, 7, 7]])
  349. """
  350. if min_distance < 0:
  351. raise ValueError(f"min_distance must be >= 0, was {min_distance}")
  352. if not np.issubdtype(label_image.dtype, np.integer):
  353. raise ValueError(
  354. f"`label_image` must be of integer dtype, got {label_image.dtype}"
  355. )
  356. if out is None:
  357. out = label_image.copy(order="C")
  358. elif out is not label_image:
  359. out[:] = label_image
  360. # May create a copy if order is not C, account for that later
  361. out_raveled = out.ravel(order="C")
  362. if spacing is not None:
  363. spacing = np.array(spacing)
  364. if spacing.shape != (out.ndim,) or spacing.min() <= 0:
  365. raise ValueError(
  366. "`spacing` must contain exactly one positive factor "
  367. "for each dimension of `label_image`"
  368. )
  369. indices = np.flatnonzero(out_raveled)
  370. # Optimization: Split indices into those on the object boundaries and inner
  371. # ones. The KDTree is built only from the boundary indices, which reduces
  372. # the size of the critical loop significantly! Remaining indices are only
  373. # used to remove the inner parts of objects as well.
  374. if (spacing is None or np.all(spacing[0] == spacing)) and p_norm <= 2:
  375. # For unity spacing we can make the borders more sparse by using a
  376. # lower connectivity
  377. footprint = ndi.generate_binary_structure(out.ndim, 1)
  378. else:
  379. footprint = ndi.generate_binary_structure(out.ndim, out.ndim)
  380. border = (
  381. ndi.maximum_filter(out, footprint=footprint)
  382. != ndi.minimum_filter(out, footprint=footprint)
  383. ).ravel()[indices]
  384. border_indices = indices[border]
  385. inner_indices = indices[~border]
  386. if border_indices.size == 0:
  387. # Image without any or only one object, return early
  388. return out
  389. # Sort by label ID first, so that IDs of the same object are contiguous
  390. # in the sorted index. This allows fast discovery of the whole object by
  391. # simple iteration up or down the index!
  392. border_indices = border_indices[np.argsort(out_raveled[border_indices])]
  393. inner_indices = inner_indices[np.argsort(out_raveled[inner_indices])]
  394. if priority is None:
  395. if not np.can_cast(out.dtype, np.intp, casting="safe"):
  396. # bincount expects intp (32-bit) on WASM or i386, so down-cast to that
  397. priority = np.bincount(out_raveled.astype(np.intp, copy=False))
  398. else:
  399. priority = np.bincount(out_raveled)
  400. # `priority` can only be indexed by positive object IDs,
  401. # `border_indices` contains all unique sorted IDs so check the lowest / first
  402. smallest_id = out_raveled[border_indices[0]]
  403. if smallest_id < 0:
  404. raise ValueError(f"found object with negative ID {smallest_id!r}")
  405. try:
  406. # Sort by priority second using a stable sort to preserve the contiguous
  407. # sorting of objects. Because each pixel in an object has the same
  408. # priority we don't need to worry about separating objects.
  409. border_indices = border_indices[
  410. np.argsort(priority[out_raveled[border_indices]], kind="stable")[::-1]
  411. ]
  412. except IndexError as error:
  413. # Use np.amax only for the exception path to provide a nicer error message
  414. expected_shape = (np.amax(out_raveled) + 1,)
  415. if priority.shape != expected_shape:
  416. raise ValueError(
  417. "shape of `priority` must be (np.amax(label_image) + 1,), "
  418. f"expected {expected_shape}, got {priority.shape} instead"
  419. ) from error
  420. else:
  421. raise
  422. # Construct kd-tree from unraveled border indices (optionally scale by `spacing`)
  423. unraveled_indices = np.unravel_index(border_indices, out.shape)
  424. if spacing is not None:
  425. unraveled_indices = tuple(
  426. unraveled_indices[dim] * spacing[dim] for dim in range(out.ndim)
  427. )
  428. kdtree = cKDTree(data=np.asarray(unraveled_indices, dtype=np.float64).T)
  429. _remove_objects_by_distance(
  430. out=out_raveled,
  431. border_indices=border_indices,
  432. inner_indices=inner_indices,
  433. kdtree=kdtree,
  434. min_distance=min_distance,
  435. p_norm=p_norm,
  436. shape=label_image.shape,
  437. )
  438. if out_raveled.base is not out:
  439. # `out_raveled` is a copy, re-assign
  440. out[:] = out_raveled.reshape(out.shape)
  441. return out