extrema.py 20 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549
  1. """extrema.py - local minima and maxima
  2. This module provides functions to find local maxima and minima of an image.
  3. Here, local maxima (minima) are defined as connected sets of pixels with equal
  4. gray level which is strictly greater (smaller) than the gray level of all
  5. pixels in direct neighborhood of the connected set. In addition, the module
  6. provides the related functions h-maxima and h-minima.
  7. Soille, P. (2003). Morphological Image Analysis: Principles and Applications
  8. (2nd ed.), Chapter 6. Springer-Verlag New York, Inc.
  9. """
  10. import numpy as np
  11. from .._shared.utils import warn
  12. from ..util import dtype_limits, invert, crop
  13. from . import grayreconstruct, _util
  14. from ._extrema_cy import _local_maxima
  15. def _add_constant_clip(image, const_value):
  16. """Add constant to the image while handling overflow issues gracefully."""
  17. min_dtype, max_dtype = dtype_limits(image, clip_negative=False)
  18. if const_value > (max_dtype - min_dtype):
  19. raise ValueError(
  20. "The added constant is not compatible" "with the image data type."
  21. )
  22. result = image + const_value
  23. result[image > max_dtype - const_value] = max_dtype
  24. return result
  25. def _subtract_constant_clip(image, const_value):
  26. """Subtract constant from image while handling underflow issues."""
  27. min_dtype, max_dtype = dtype_limits(image, clip_negative=False)
  28. if const_value > (max_dtype - min_dtype):
  29. raise ValueError(
  30. "The subtracted constant is not compatible" "with the image data type."
  31. )
  32. result = image - const_value
  33. result[image < (const_value + min_dtype)] = min_dtype
  34. return result
  35. def h_maxima(image, h, footprint=None):
  36. """Determine all maxima of the image with height >= h.
  37. The local maxima are defined as connected sets of pixels with equal
  38. gray level strictly greater than the gray level of all pixels in direct
  39. neighborhood of the set.
  40. A local maximum M of height h is a local maximum for which
  41. there is at least one path joining M with an equal or higher local maximum
  42. on which the minimal value is f(M) - h (i.e. the values along the path
  43. are not decreasing by more than h with respect to the maximum's value)
  44. and no path to an equal or higher local maximum for which the minimal
  45. value is greater.
  46. The global maxima of the image are also found by this function.
  47. Parameters
  48. ----------
  49. image : ndarray
  50. The input image for which the maxima are to be calculated.
  51. h : unsigned integer
  52. The minimal height of all extracted maxima.
  53. footprint : ndarray, optional
  54. The neighborhood expressed as an n-D array of 1's and 0's.
  55. Default is the ball of radius 1 according to the maximum norm
  56. (i.e. a 3x3 square for 2D images, a 3x3x3 cube for 3D images, etc.)
  57. Returns
  58. -------
  59. h_max : ndarray
  60. The local maxima of height >= h and the global maxima.
  61. The resulting image is a binary image, where pixels belonging to
  62. the determined maxima take value 1, the others take value 0.
  63. See Also
  64. --------
  65. skimage.morphology.h_minima
  66. skimage.morphology.local_maxima
  67. skimage.morphology.local_minima
  68. References
  69. ----------
  70. .. [1] Soille, P., "Morphological Image Analysis: Principles and
  71. Applications" (Chapter 6), 2nd edition (2003), ISBN 3540429883.
  72. Examples
  73. --------
  74. >>> import numpy as np
  75. >>> from skimage.morphology import extrema
  76. We create an image (quadratic function with a maximum in the center and
  77. 4 additional constant maxima.
  78. The heights of the maxima are: 1, 21, 41, 61, 81
  79. >>> w = 10
  80. >>> x, y = np.mgrid[0:w,0:w]
  81. >>> f = 20 - 0.2*((x - w/2)**2 + (y-w/2)**2)
  82. >>> f[2:4,2:4] = 40; f[2:4,7:9] = 60; f[7:9,2:4] = 80; f[7:9,7:9] = 100
  83. >>> f = f.astype(int)
  84. We can calculate all maxima with a height of at least 40:
  85. >>> maxima = extrema.h_maxima(f, 40)
  86. The resulting image will contain 3 local maxima.
  87. """
  88. # Check for h value that is larger then range of the image. If this
  89. # is True then there are no h-maxima in the image.
  90. if h > np.ptp(image):
  91. return np.zeros(image.shape, dtype=np.uint8)
  92. # Check for floating point h value. For this to work properly
  93. # we need to explicitly convert image to float64.
  94. #
  95. # FIXME: This could give incorrect results if image is int64 and
  96. # has a very high dynamic range. The dtype of image is
  97. # changed to float64, and different integer values could
  98. # become the same float due to rounding.
  99. #
  100. # >>> ii64 = np.iinfo(np.int64)
  101. # >>> a = np.array([ii64.max, ii64.max - 2])
  102. # >>> a[0] == a[1]
  103. # False
  104. # >>> b = a.astype(np.float64)
  105. # >>> b[0] == b[1]
  106. # True
  107. #
  108. if np.issubdtype(type(h), np.floating) and np.issubdtype(image.dtype, np.integer):
  109. if (h % 1) != 0:
  110. warn(
  111. 'possible precision loss converting image to '
  112. 'floating point. To silence this warning, '
  113. 'ensure image and h have same data type.',
  114. stacklevel=2,
  115. )
  116. image = image.astype(float)
  117. else:
  118. h = image.dtype.type(h)
  119. if h == 0:
  120. raise ValueError("h = 0 is ambiguous, use local_maxima() " "instead?")
  121. if np.issubdtype(image.dtype, np.floating):
  122. # The purpose of the resolution variable is to allow for the
  123. # small rounding errors that inevitably occur when doing
  124. # floating point arithmetic. We want shifted_img to be
  125. # guaranteed to be h less than image. If we only subtract h
  126. # there may be pixels were shifted_img ends up being
  127. # slightly greater than image - h.
  128. #
  129. # The resolution is scaled based on the pixel values in the
  130. # image because floating point precision is relative. A
  131. # very large value of 1.0e10 will have a large precision,
  132. # say +-1.0e4, and a very small value of 1.0e-10 will have
  133. # a very small precision, say +-1.0e-16.
  134. #
  135. resolution = 2 * np.finfo(image.dtype).resolution * np.abs(image)
  136. shifted_img = image - h - resolution
  137. else:
  138. shifted_img = _subtract_constant_clip(image, h)
  139. rec_img = grayreconstruct.reconstruction(
  140. shifted_img, image, method='dilation', footprint=footprint
  141. )
  142. residue_img = image - rec_img
  143. return (residue_img >= h).astype(np.uint8)
  144. def h_minima(image, h, footprint=None):
  145. """Determine all minima of the image with depth >= h.
  146. The local minima are defined as connected sets of pixels with equal
  147. gray level strictly smaller than the gray levels of all pixels in direct
  148. neighborhood of the set.
  149. A local minimum M of depth h is a local minimum for which
  150. there is at least one path joining M with an equal or lower local minimum
  151. on which the maximal value is f(M) + h (i.e. the values along the path
  152. are not increasing by more than h with respect to the minimum's value)
  153. and no path to an equal or lower local minimum for which the maximal
  154. value is smaller.
  155. The global minima of the image are also found by this function.
  156. Parameters
  157. ----------
  158. image : ndarray
  159. The input image for which the minima are to be calculated.
  160. h : unsigned integer
  161. The minimal depth of all extracted minima.
  162. footprint : ndarray, optional
  163. The neighborhood expressed as an n-D array of 1's and 0's.
  164. Default is the ball of radius 1 according to the maximum norm
  165. (i.e. a 3x3 square for 2D images, a 3x3x3 cube for 3D images, etc.)
  166. Returns
  167. -------
  168. h_min : ndarray
  169. The local minima of depth >= h and the global minima.
  170. The resulting image is a binary image, where pixels belonging to
  171. the determined minima take value 1, the others take value 0.
  172. See Also
  173. --------
  174. skimage.morphology.h_maxima
  175. skimage.morphology.local_maxima
  176. skimage.morphology.local_minima
  177. References
  178. ----------
  179. .. [1] Soille, P., "Morphological Image Analysis: Principles and
  180. Applications" (Chapter 6), 2nd edition (2003), ISBN 3540429883.
  181. Examples
  182. --------
  183. >>> import numpy as np
  184. >>> from skimage.morphology import extrema
  185. We create an image (quadratic function with a minimum in the center and
  186. 4 additional constant maxima.
  187. The depth of the minima are: 1, 21, 41, 61, 81
  188. >>> w = 10
  189. >>> x, y = np.mgrid[0:w,0:w]
  190. >>> f = 180 + 0.2*((x - w/2)**2 + (y-w/2)**2)
  191. >>> f[2:4,2:4] = 160; f[2:4,7:9] = 140; f[7:9,2:4] = 120; f[7:9,7:9] = 100
  192. >>> f = f.astype(int)
  193. We can calculate all minima with a depth of at least 40:
  194. >>> minima = extrema.h_minima(f, 40)
  195. The resulting image will contain 3 local minima.
  196. """
  197. if h > np.ptp(image):
  198. return np.zeros(image.shape, dtype=np.uint8)
  199. if np.issubdtype(type(h), np.floating) and np.issubdtype(image.dtype, np.integer):
  200. if (h % 1) != 0:
  201. warn(
  202. 'possible precision loss converting image to '
  203. 'floating point. To silence this warning, '
  204. 'ensure image and h have same data type.',
  205. stacklevel=2,
  206. )
  207. image = image.astype(float)
  208. else:
  209. h = image.dtype.type(h)
  210. if h == 0:
  211. raise ValueError("h = 0 is ambiguous, use local_minima() " "instead?")
  212. if np.issubdtype(image.dtype, np.floating):
  213. resolution = 2 * np.finfo(image.dtype).resolution * np.abs(image)
  214. shifted_img = image + h + resolution
  215. else:
  216. shifted_img = _add_constant_clip(image, h)
  217. rec_img = grayreconstruct.reconstruction(
  218. shifted_img, image, method='erosion', footprint=footprint
  219. )
  220. residue_img = rec_img - image
  221. return (residue_img >= h).astype(np.uint8)
  222. def local_maxima(
  223. image, footprint=None, connectivity=None, indices=False, allow_borders=True
  224. ):
  225. """Find local maxima of n-dimensional array.
  226. The local maxima are defined as connected sets of pixels with equal gray
  227. level (plateaus) strictly greater than the gray levels of all pixels in the
  228. neighborhood.
  229. Parameters
  230. ----------
  231. image : ndarray
  232. An n-dimensional array.
  233. footprint : ndarray, optional
  234. The footprint (structuring element) used to determine the neighborhood
  235. of each evaluated pixel (``True`` denotes a connected pixel). It must
  236. be a boolean array and have the same number of dimensions as `image`.
  237. If neither `footprint` nor `connectivity` are given, all adjacent
  238. pixels are considered as part of the neighborhood.
  239. connectivity : int, optional
  240. A number used to determine the neighborhood of each evaluated pixel.
  241. Adjacent pixels whose squared distance from the center is less than or
  242. equal to `connectivity` are considered neighbors. Ignored if
  243. `footprint` is not None.
  244. indices : bool, optional
  245. If True, the output will be a tuple of one-dimensional arrays
  246. representing the indices of local maxima in each dimension. If False,
  247. the output will be a boolean array with the same shape as `image`.
  248. allow_borders : bool, optional
  249. If true, plateaus that touch the image border are valid maxima.
  250. Returns
  251. -------
  252. maxima : ndarray or tuple[ndarray]
  253. If `indices` is false, a boolean array with the same shape as `image`
  254. is returned with ``True`` indicating the position of local maxima
  255. (``False`` otherwise). If `indices` is true, a tuple of one-dimensional
  256. arrays containing the coordinates (indices) of all found maxima.
  257. Warns
  258. -----
  259. UserWarning
  260. If `allow_borders` is false and any dimension of the given `image` is
  261. shorter than 3 samples, maxima can't exist and a warning is shown.
  262. See Also
  263. --------
  264. skimage.morphology.local_minima
  265. skimage.morphology.h_maxima
  266. skimage.morphology.h_minima
  267. Notes
  268. -----
  269. This function operates on the following ideas:
  270. 1. Make a first pass over the image's last dimension and flag candidates
  271. for local maxima by comparing pixels in only one direction.
  272. If the pixels aren't connected in the last dimension all pixels are
  273. flagged as candidates instead.
  274. For each candidate:
  275. 2. Perform a flood-fill to find all connected pixels that have the same
  276. gray value and are part of the plateau.
  277. 3. Consider the connected neighborhood of a plateau: if no bordering sample
  278. has a higher gray level, mark the plateau as a definite local maximum.
  279. Examples
  280. --------
  281. >>> from skimage.morphology import local_maxima
  282. >>> image = np.zeros((4, 7), dtype=int)
  283. >>> image[1:3, 1:3] = 1
  284. >>> image[3, 0] = 1
  285. >>> image[1:3, 4:6] = 2
  286. >>> image[3, 6] = 3
  287. >>> image
  288. array([[0, 0, 0, 0, 0, 0, 0],
  289. [0, 1, 1, 0, 2, 2, 0],
  290. [0, 1, 1, 0, 2, 2, 0],
  291. [1, 0, 0, 0, 0, 0, 3]])
  292. Find local maxima by comparing to all neighboring pixels (maximal
  293. connectivity):
  294. >>> local_maxima(image)
  295. array([[False, False, False, False, False, False, False],
  296. [False, True, True, False, False, False, False],
  297. [False, True, True, False, False, False, False],
  298. [ True, False, False, False, False, False, True]])
  299. >>> local_maxima(image, indices=True)
  300. (array([1, 1, 2, 2, 3, 3]), array([1, 2, 1, 2, 0, 6]))
  301. Find local maxima without comparing to diagonal pixels (connectivity 1):
  302. >>> local_maxima(image, connectivity=1)
  303. array([[False, False, False, False, False, False, False],
  304. [False, True, True, False, True, True, False],
  305. [False, True, True, False, True, True, False],
  306. [ True, False, False, False, False, False, True]])
  307. and exclude maxima that border the image edge:
  308. >>> local_maxima(image, connectivity=1, allow_borders=False)
  309. array([[False, False, False, False, False, False, False],
  310. [False, True, True, False, True, True, False],
  311. [False, True, True, False, True, True, False],
  312. [False, False, False, False, False, False, False]])
  313. """
  314. image = np.asarray(image, order="C")
  315. if image.size == 0:
  316. # Return early for empty input
  317. if indices:
  318. # Make sure that output is a tuple of 1 empty array per dimension
  319. return np.nonzero(image)
  320. else:
  321. return np.zeros(image.shape, dtype=bool)
  322. if allow_borders:
  323. # Ensure that local maxima are always at least one smaller sample away
  324. # from the image border
  325. image = np.pad(image, 1, mode='constant', constant_values=image.min())
  326. # Array of flags used to store the state of each pixel during evaluation.
  327. # See _extrema_cy.pyx for their meaning
  328. flags = np.zeros(image.shape, dtype=np.uint8)
  329. _util._set_border_values(flags, value=3)
  330. if any(s < 3 for s in image.shape):
  331. # Warn and skip if any dimension is smaller than 3
  332. # -> no maxima can exist & footprint can't be applied
  333. warn(
  334. "maxima can't exist for an image with any dimension smaller 3 "
  335. "if borders aren't allowed",
  336. stacklevel=3,
  337. )
  338. else:
  339. footprint = _util._resolve_neighborhood(footprint, connectivity, image.ndim)
  340. neighbor_offsets = _util._offsets_to_raveled_neighbors(
  341. image.shape, footprint, center=((1,) * image.ndim)
  342. )
  343. try:
  344. _local_maxima(image.ravel(), flags.ravel(), neighbor_offsets)
  345. except TypeError:
  346. if image.dtype == np.float16:
  347. # Provide the user with clearer error message
  348. raise TypeError(
  349. "dtype of `image` is float16 which is not "
  350. "supported, try upcasting to float32"
  351. )
  352. else:
  353. raise # Otherwise raise original message
  354. if allow_borders:
  355. # Revert padding performed at the beginning of the function
  356. flags = crop(flags, 1)
  357. else:
  358. # No padding was performed but set edge values back to 0
  359. _util._set_border_values(flags, value=0)
  360. if indices:
  361. return np.nonzero(flags)
  362. else:
  363. return flags.view(bool)
  364. def local_minima(
  365. image, footprint=None, connectivity=None, indices=False, allow_borders=True
  366. ):
  367. """Find local minima of n-dimensional array.
  368. The local minima are defined as connected sets of pixels with equal gray
  369. level (plateaus) strictly smaller than the gray levels of all pixels in the
  370. neighborhood.
  371. Parameters
  372. ----------
  373. image : ndarray
  374. An n-dimensional array.
  375. footprint : ndarray, optional
  376. The footprint (structuring element) used to determine the neighborhood
  377. of each evaluated pixel (``True`` denotes a connected pixel). It must
  378. be a boolean array and have the same number of dimensions as `image`.
  379. If neither `footprint` nor `connectivity` are given, all adjacent
  380. pixels are considered as part of the neighborhood.
  381. connectivity : int, optional
  382. A number used to determine the neighborhood of each evaluated pixel.
  383. Adjacent pixels whose squared distance from the center is less than or
  384. equal to `connectivity` are considered neighbors. Ignored if
  385. `footprint` is not None.
  386. indices : bool, optional
  387. If True, the output will be a tuple of one-dimensional arrays
  388. representing the indices of local minima in each dimension. If False,
  389. the output will be a boolean array with the same shape as `image`.
  390. allow_borders : bool, optional
  391. If true, plateaus that touch the image border are valid minima.
  392. Returns
  393. -------
  394. minima : ndarray or tuple[ndarray]
  395. If `indices` is false, a boolean array with the same shape as `image`
  396. is returned with ``True`` indicating the position of local minima
  397. (``False`` otherwise). If `indices` is true, a tuple of one-dimensional
  398. arrays containing the coordinates (indices) of all found minima.
  399. See Also
  400. --------
  401. skimage.morphology.local_maxima
  402. skimage.morphology.h_maxima
  403. skimage.morphology.h_minima
  404. Notes
  405. -----
  406. This function operates on the following ideas:
  407. 1. Make a first pass over the image's last dimension and flag candidates
  408. for local minima by comparing pixels in only one direction.
  409. If the pixels aren't connected in the last dimension all pixels are
  410. flagged as candidates instead.
  411. For each candidate:
  412. 2. Perform a flood-fill to find all connected pixels that have the same
  413. gray value and are part of the plateau.
  414. 3. Consider the connected neighborhood of a plateau: if no bordering sample
  415. has a smaller gray level, mark the plateau as a definite local minimum.
  416. Examples
  417. --------
  418. >>> from skimage.morphology import local_minima
  419. >>> image = np.zeros((4, 7), dtype=int)
  420. >>> image[1:3, 1:3] = -1
  421. >>> image[3, 0] = -1
  422. >>> image[1:3, 4:6] = -2
  423. >>> image[3, 6] = -3
  424. >>> image
  425. array([[ 0, 0, 0, 0, 0, 0, 0],
  426. [ 0, -1, -1, 0, -2, -2, 0],
  427. [ 0, -1, -1, 0, -2, -2, 0],
  428. [-1, 0, 0, 0, 0, 0, -3]])
  429. Find local minima by comparing to all neighboring pixels (maximal
  430. connectivity):
  431. >>> local_minima(image)
  432. array([[False, False, False, False, False, False, False],
  433. [False, True, True, False, False, False, False],
  434. [False, True, True, False, False, False, False],
  435. [ True, False, False, False, False, False, True]])
  436. >>> local_minima(image, indices=True)
  437. (array([1, 1, 2, 2, 3, 3]), array([1, 2, 1, 2, 0, 6]))
  438. Find local minima without comparing to diagonal pixels (connectivity 1):
  439. >>> local_minima(image, connectivity=1)
  440. array([[False, False, False, False, False, False, False],
  441. [False, True, True, False, True, True, False],
  442. [False, True, True, False, True, True, False],
  443. [ True, False, False, False, False, False, True]])
  444. and exclude minima that border the image edge:
  445. >>> local_minima(image, connectivity=1, allow_borders=False)
  446. array([[False, False, False, False, False, False, False],
  447. [False, True, True, False, True, True, False],
  448. [False, True, True, False, True, True, False],
  449. [False, False, False, False, False, False, False]])
  450. """
  451. return local_maxima(
  452. image=invert(image, signed_float=True),
  453. footprint=footprint,
  454. connectivity=connectivity,
  455. indices=indices,
  456. allow_borders=allow_borders,
  457. )