max_tree.py 26 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700
  1. """max_tree.py - max_tree representation of images.
  2. This module provides operators based on the max-tree representation of images.
  3. A grayscale image can be seen as a pile of nested sets, each of which is the
  4. result of a threshold operation. These sets can be efficiently represented by
  5. max-trees, where the inclusion relation between connected components at
  6. different levels are represented by parent-child relationships.
  7. These representations allow efficient implementations of many algorithms, such
  8. as attribute operators. Unlike morphological openings and closings, these
  9. operators do not require a fixed footprint, but rather act with a flexible
  10. footprint that meets a certain criterion.
  11. This implementation provides functions for:
  12. 1. max-tree generation
  13. 2. area openings / closings
  14. 3. diameter openings / closings
  15. 4. local maxima
  16. References:
  17. .. [1] Salembier, P., Oliveras, A., & Garrido, L. (1998). Antiextensive
  18. Connected Operators for Image and Sequence Processing.
  19. IEEE Transactions on Image Processing, 7(4), 555-570.
  20. :DOI:`10.1109/83.663500`
  21. .. [2] Berger, C., Geraud, T., Levillain, R., Widynski, N., Baillard, A.,
  22. Bertin, E. (2007). Effective Component Tree Computation with
  23. Application to Pattern Recognition in Astronomical Imaging.
  24. In International Conference on Image Processing (ICIP) (pp. 41-44).
  25. :DOI:`10.1109/ICIP.2007.4379949`
  26. .. [3] Najman, L., & Couprie, M. (2006). Building the component tree in
  27. quasi-linear time. IEEE Transactions on Image Processing, 15(11),
  28. 3531-3539.
  29. :DOI:`10.1109/TIP.2006.877518`
  30. .. [4] Carlinet, E., & Geraud, T. (2014). A Comparative Review of
  31. Component Tree Computation Algorithms. IEEE Transactions on Image
  32. Processing, 23(9), 3885-3895.
  33. :DOI:`10.1109/TIP.2014.2336551`
  34. """
  35. import numpy as np
  36. from ._util import _validate_connectivity, _offsets_to_raveled_neighbors
  37. from ..util import invert
  38. from . import _max_tree
  39. unsigned_int_types = [np.uint8, np.uint16, np.uint32, np.uint64]
  40. signed_int_types = [np.int8, np.int16, np.int32, np.int64]
  41. signed_float_types = [np.float16, np.float32, np.float64]
  42. # building the max tree.
  43. def max_tree(image, connectivity=1):
  44. """Build the max tree from an image.
  45. Component trees represent the hierarchical structure of the connected
  46. components resulting from sequential thresholding operations applied to an
  47. image. A connected component at one level is parent of a component at a
  48. higher level if the latter is included in the first. A max-tree is an
  49. efficient representation of a component tree. A connected component at
  50. one level is represented by one reference pixel at this level, which is
  51. parent to all other pixels at that level and to the reference pixel at the
  52. level above. The max-tree is the basis for many morphological operators,
  53. namely connected operators.
  54. Parameters
  55. ----------
  56. image : ndarray
  57. The input image for which the max-tree is to be calculated.
  58. This image can be of any type.
  59. connectivity : unsigned int, optional
  60. The neighborhood connectivity. The integer represents the maximum
  61. number of orthogonal steps to reach a neighbor. In 2D, it is 1 for
  62. a 4-neighborhood and 2 for a 8-neighborhood. Default value is 1.
  63. Returns
  64. -------
  65. parent : ndarray, int64
  66. Array of same shape as image. The value of each pixel is the index of
  67. its parent in the ravelled array.
  68. tree_traverser : 1D array, int64
  69. The ordered pixel indices (referring to the ravelled array). The pixels
  70. are ordered such that every pixel is preceded by its parent (except for
  71. the root which has no parent).
  72. References
  73. ----------
  74. .. [1] Salembier, P., Oliveras, A., & Garrido, L. (1998). Antiextensive
  75. Connected Operators for Image and Sequence Processing.
  76. IEEE Transactions on Image Processing, 7(4), 555-570.
  77. :DOI:`10.1109/83.663500`
  78. .. [2] Berger, C., Geraud, T., Levillain, R., Widynski, N., Baillard, A.,
  79. Bertin, E. (2007). Effective Component Tree Computation with
  80. Application to Pattern Recognition in Astronomical Imaging.
  81. In International Conference on Image Processing (ICIP) (pp. 41-44).
  82. :DOI:`10.1109/ICIP.2007.4379949`
  83. .. [3] Najman, L., & Couprie, M. (2006). Building the component tree in
  84. quasi-linear time. IEEE Transactions on Image Processing, 15(11),
  85. 3531-3539.
  86. :DOI:`10.1109/TIP.2006.877518`
  87. .. [4] Carlinet, E., & Geraud, T. (2014). A Comparative Review of
  88. Component Tree Computation Algorithms. IEEE Transactions on Image
  89. Processing, 23(9), 3885-3895.
  90. :DOI:`10.1109/TIP.2014.2336551`
  91. Examples
  92. --------
  93. We create a small sample image (Figure 1 from [4]) and build the max-tree.
  94. >>> image = np.array([[15, 13, 16], [12, 12, 10], [16, 12, 14]])
  95. >>> P, S = max_tree(image, connectivity=2)
  96. """
  97. # User defined masks are not allowed, as there might be more than one
  98. # connected component in the mask (and therefore not a single tree that
  99. # represents the image). Mask here is an image that is 0 on the border
  100. # and 1 everywhere else.
  101. mask = np.ones(image.shape)
  102. for k in range(len(image.shape)):
  103. np.moveaxis(mask, k, 0)[0] = 0
  104. np.moveaxis(mask, k, 0)[-1] = 0
  105. neighbors, offset = _validate_connectivity(image.ndim, connectivity, offset=None)
  106. # initialization of the parent image
  107. parent = np.zeros(image.shape, dtype=np.int64)
  108. # flat_neighborhood contains a list of offsets allowing one to find the
  109. # neighbors in the ravelled image.
  110. flat_neighborhood = _offsets_to_raveled_neighbors(
  111. image.shape, neighbors, offset
  112. ).astype(np.int32)
  113. # pixels need to be sorted according to their gray level.
  114. tree_traverser = np.argsort(image.ravel(), kind="stable").astype(np.int64)
  115. # call of cython function.
  116. _max_tree._max_tree(
  117. image.ravel(),
  118. mask.ravel().astype(np.uint8),
  119. flat_neighborhood,
  120. offset.astype(np.int32),
  121. np.array(image.shape, dtype=np.int32),
  122. parent.ravel(),
  123. tree_traverser,
  124. )
  125. return parent, tree_traverser
  126. def area_opening(
  127. image, area_threshold=64, connectivity=1, parent=None, tree_traverser=None
  128. ):
  129. """Perform an area opening of the image.
  130. Area opening removes all bright structures of an image with
  131. a surface smaller than area_threshold.
  132. The output image is thus the largest image smaller than the input
  133. for which all local maxima have at least a surface of
  134. area_threshold pixels.
  135. Area openings are similar to morphological openings, but
  136. they do not use a fixed footprint, but rather a deformable
  137. one, with surface = area_threshold. Consequently, the area_opening
  138. with area_threshold=1 is the identity.
  139. In the binary case, area openings are equivalent to
  140. remove_small_objects; this operator is thus extended to gray-level images.
  141. Technically, this operator is based on the max-tree representation of
  142. the image.
  143. Parameters
  144. ----------
  145. image : ndarray
  146. The input image for which the area_opening is to be calculated.
  147. This image can be of any type.
  148. area_threshold : unsigned int
  149. The size parameter (number of pixels). The default value is arbitrarily
  150. chosen to be 64.
  151. connectivity : unsigned int, optional
  152. The neighborhood connectivity. The integer represents the maximum
  153. number of orthogonal steps to reach a neighbor. In 2D, it is 1 for
  154. a 4-neighborhood and 2 for a 8-neighborhood. Default value is 1.
  155. parent : ndarray, int64, optional
  156. Parent image representing the max tree of the image. The
  157. value of each pixel is the index of its parent in the ravelled array.
  158. tree_traverser : 1D array, int64, optional
  159. The ordered pixel indices (referring to the ravelled array). The pixels
  160. are ordered such that every pixel is preceded by its parent (except for
  161. the root which has no parent).
  162. Returns
  163. -------
  164. output : ndarray
  165. Output image of the same shape and type as the input image.
  166. See Also
  167. --------
  168. skimage.morphology.area_closing
  169. skimage.morphology.diameter_opening
  170. skimage.morphology.diameter_closing
  171. skimage.morphology.max_tree
  172. skimage.morphology.remove_small_objects
  173. skimage.morphology.remove_small_holes
  174. References
  175. ----------
  176. .. [1] Vincent L., Proc. "Grayscale area openings and closings,
  177. their efficient implementation and applications",
  178. EURASIP Workshop on Mathematical Morphology and its
  179. Applications to Signal Processing, Barcelona, Spain, pp.22-27,
  180. May 1993.
  181. .. [2] Soille, P., "Morphological Image Analysis: Principles and
  182. Applications" (Chapter 6), 2nd edition (2003), ISBN 3540429883.
  183. :DOI:`10.1007/978-3-662-05088-0`
  184. .. [3] Salembier, P., Oliveras, A., & Garrido, L. (1998). Antiextensive
  185. Connected Operators for Image and Sequence Processing.
  186. IEEE Transactions on Image Processing, 7(4), 555-570.
  187. :DOI:`10.1109/83.663500`
  188. .. [4] Najman, L., & Couprie, M. (2006). Building the component tree in
  189. quasi-linear time. IEEE Transactions on Image Processing, 15(11),
  190. 3531-3539.
  191. :DOI:`10.1109/TIP.2006.877518`
  192. .. [5] Carlinet, E., & Geraud, T. (2014). A Comparative Review of
  193. Component Tree Computation Algorithms. IEEE Transactions on Image
  194. Processing, 23(9), 3885-3895.
  195. :DOI:`10.1109/TIP.2014.2336551`
  196. Examples
  197. --------
  198. We create an image (quadratic function with a maximum in the center and
  199. 4 additional local maxima.
  200. >>> w = 12
  201. >>> x, y = np.mgrid[0:w,0:w]
  202. >>> f = 20 - 0.2*((x - w/2)**2 + (y-w/2)**2)
  203. >>> f[2:3,1:5] = 40; f[2:4,9:11] = 60; f[9:11,2:4] = 80
  204. >>> f[9:10,9:11] = 100; f[10,10] = 100
  205. >>> f = f.astype(int)
  206. We can calculate the area opening:
  207. >>> open = area_opening(f, 8, connectivity=1)
  208. The peaks with a surface smaller than 8 are removed.
  209. """
  210. output = image.copy()
  211. if parent is None or tree_traverser is None:
  212. parent, tree_traverser = max_tree(image, connectivity)
  213. area = _max_tree._compute_area(image.ravel(), parent.ravel(), tree_traverser)
  214. _max_tree._direct_filter(
  215. image.ravel(),
  216. output.ravel(),
  217. parent.ravel(),
  218. tree_traverser,
  219. area,
  220. area_threshold,
  221. )
  222. return output
  223. def diameter_opening(
  224. image, diameter_threshold=8, connectivity=1, parent=None, tree_traverser=None
  225. ):
  226. """Perform a diameter opening of the image.
  227. Diameter opening removes all bright structures of an image with
  228. maximal extension smaller than diameter_threshold. The maximal
  229. extension is defined as the maximal extension of the bounding box.
  230. The operator is also called Bounding Box Opening. In practice,
  231. the result is similar to a morphological opening, but long and thin
  232. structures are not removed.
  233. Technically, this operator is based on the max-tree representation of
  234. the image.
  235. Parameters
  236. ----------
  237. image : ndarray
  238. The input image for which the area_opening is to be calculated.
  239. This image can be of any type.
  240. diameter_threshold : unsigned int
  241. The maximal extension parameter (number of pixels). The default value
  242. is 8.
  243. connectivity : unsigned int, optional
  244. The neighborhood connectivity. The integer represents the maximum
  245. number of orthogonal steps to reach a neighbor. In 2D, it is 1 for
  246. a 4-neighborhood and 2 for a 8-neighborhood. Default value is 1.
  247. parent : ndarray, int64, optional
  248. Parent image representing the max tree of the image. The
  249. value of each pixel is the index of its parent in the ravelled array.
  250. tree_traverser : 1D array, int64, optional
  251. The ordered pixel indices (referring to the ravelled array). The pixels
  252. are ordered such that every pixel is preceded by its parent (except for
  253. the root which has no parent).
  254. Returns
  255. -------
  256. output : ndarray
  257. Output image of the same shape and type as the input image.
  258. See Also
  259. --------
  260. skimage.morphology.area_opening
  261. skimage.morphology.area_closing
  262. skimage.morphology.diameter_closing
  263. skimage.morphology.max_tree
  264. References
  265. ----------
  266. .. [1] Walter, T., & Klein, J.-C. (2002). Automatic Detection of
  267. Microaneurysms in Color Fundus Images of the Human Retina by Means
  268. of the Bounding Box Closing. In A. Colosimo, P. Sirabella,
  269. A. Giuliani (Eds.), Medical Data Analysis. Lecture Notes in Computer
  270. Science, vol 2526, pp. 210-220. Springer Berlin Heidelberg.
  271. :DOI:`10.1007/3-540-36104-9_23`
  272. .. [2] Carlinet, E., & Geraud, T. (2014). A Comparative Review of
  273. Component Tree Computation Algorithms. IEEE Transactions on Image
  274. Processing, 23(9), 3885-3895.
  275. :DOI:`10.1109/TIP.2014.2336551`
  276. Examples
  277. --------
  278. We create an image (quadratic function with a maximum in the center and
  279. 4 additional local maxima.
  280. >>> w = 12
  281. >>> x, y = np.mgrid[0:w,0:w]
  282. >>> f = 20 - 0.2*((x - w/2)**2 + (y-w/2)**2)
  283. >>> f[2:3,1:5] = 40; f[2:4,9:11] = 60; f[9:11,2:4] = 80
  284. >>> f[9:10,9:11] = 100; f[10,10] = 100
  285. >>> f = f.astype(int)
  286. We can calculate the diameter opening:
  287. >>> open = diameter_opening(f, 3, connectivity=1)
  288. The peaks with a maximal extension of 2 or less are removed.
  289. The remaining peaks have all a maximal extension of at least 3.
  290. """
  291. output = image.copy()
  292. if parent is None or tree_traverser is None:
  293. parent, tree_traverser = max_tree(image, connectivity)
  294. diam = _max_tree._compute_extension(
  295. image.ravel(),
  296. np.array(image.shape, dtype=np.int32),
  297. parent.ravel(),
  298. tree_traverser,
  299. )
  300. _max_tree._direct_filter(
  301. image.ravel(),
  302. output.ravel(),
  303. parent.ravel(),
  304. tree_traverser,
  305. diam,
  306. diameter_threshold,
  307. )
  308. return output
  309. def area_closing(
  310. image, area_threshold=64, connectivity=1, parent=None, tree_traverser=None
  311. ):
  312. """Perform an area closing of the image.
  313. Area closing removes all dark structures of an image with
  314. a surface smaller than area_threshold.
  315. The output image is larger than or equal to the input image
  316. for every pixel and all local minima have at least a surface of
  317. area_threshold pixels.
  318. Area closings are similar to morphological closings, but
  319. they do not use a fixed footprint, but rather a deformable
  320. one, with surface = area_threshold.
  321. In the binary case, area closings are equivalent to
  322. remove_small_holes; this operator is thus extended to gray-level images.
  323. Technically, this operator is based on the max-tree representation of
  324. the image.
  325. Parameters
  326. ----------
  327. image : ndarray
  328. The input image for which the area_closing is to be calculated.
  329. This image can be of any type.
  330. area_threshold : unsigned int
  331. The size parameter (number of pixels). The default value is arbitrarily
  332. chosen to be 64.
  333. connectivity : unsigned int, optional
  334. The neighborhood connectivity. The integer represents the maximum
  335. number of orthogonal steps to reach a neighbor. In 2D, it is 1 for
  336. a 4-neighborhood and 2 for a 8-neighborhood. Default value is 1.
  337. parent : ndarray, int64, optional
  338. Parent image representing the max tree of the inverted image. The
  339. value of each pixel is the index of its parent in the ravelled array.
  340. See Note for further details.
  341. tree_traverser : 1D array, int64, optional
  342. The ordered pixel indices (referring to the ravelled array). The pixels
  343. are ordered such that every pixel is preceded by its parent (except for
  344. the root which has no parent).
  345. Returns
  346. -------
  347. output : ndarray
  348. Output image of the same shape and type as input image.
  349. See Also
  350. --------
  351. skimage.morphology.area_opening
  352. skimage.morphology.diameter_opening
  353. skimage.morphology.diameter_closing
  354. skimage.morphology.max_tree
  355. skimage.morphology.remove_small_objects
  356. skimage.morphology.remove_small_holes
  357. References
  358. ----------
  359. .. [1] Vincent L., Proc. "Grayscale area openings and closings,
  360. their efficient implementation and applications",
  361. EURASIP Workshop on Mathematical Morphology and its
  362. Applications to Signal Processing, Barcelona, Spain, pp.22-27,
  363. May 1993.
  364. .. [2] Soille, P., "Morphological Image Analysis: Principles and
  365. Applications" (Chapter 6), 2nd edition (2003), ISBN 3540429883.
  366. :DOI:`10.1007/978-3-662-05088-0`
  367. .. [3] Salembier, P., Oliveras, A., & Garrido, L. (1998). Antiextensive
  368. Connected Operators for Image and Sequence Processing.
  369. IEEE Transactions on Image Processing, 7(4), 555-570.
  370. :DOI:`10.1109/83.663500`
  371. .. [4] Najman, L., & Couprie, M. (2006). Building the component tree in
  372. quasi-linear time. IEEE Transactions on Image Processing, 15(11),
  373. 3531-3539.
  374. :DOI:`10.1109/TIP.2006.877518`
  375. .. [5] Carlinet, E., & Geraud, T. (2014). A Comparative Review of
  376. Component Tree Computation Algorithms. IEEE Transactions on Image
  377. Processing, 23(9), 3885-3895.
  378. :DOI:`10.1109/TIP.2014.2336551`
  379. Examples
  380. --------
  381. We create an image (quadratic function with a minimum in the center and
  382. 4 additional local minima.
  383. >>> w = 12
  384. >>> x, y = np.mgrid[0:w,0:w]
  385. >>> f = 180 + 0.2*((x - w/2)**2 + (y-w/2)**2)
  386. >>> f[2:3,1:5] = 160; f[2:4,9:11] = 140; f[9:11,2:4] = 120
  387. >>> f[9:10,9:11] = 100; f[10,10] = 100
  388. >>> f = f.astype(int)
  389. We can calculate the area closing:
  390. >>> closed = area_closing(f, 8, connectivity=1)
  391. All small minima are removed, and the remaining minima have at least
  392. a size of 8.
  393. Notes
  394. -----
  395. If a max-tree representation (parent and tree_traverser) are given to the
  396. function, they must be calculated from the inverted image for this
  397. function, i.e.:
  398. >>> P, S = max_tree(invert(f))
  399. >>> closed = diameter_closing(f, 3, parent=P, tree_traverser=S)
  400. """
  401. # inversion of the input image
  402. image_inv = invert(image)
  403. output = image_inv.copy()
  404. if parent is None or tree_traverser is None:
  405. parent, tree_traverser = max_tree(image_inv, connectivity)
  406. area = _max_tree._compute_area(image_inv.ravel(), parent.ravel(), tree_traverser)
  407. _max_tree._direct_filter(
  408. image_inv.ravel(),
  409. output.ravel(),
  410. parent.ravel(),
  411. tree_traverser,
  412. area,
  413. area_threshold,
  414. )
  415. # inversion of the output image
  416. output = invert(output)
  417. return output
  418. def diameter_closing(
  419. image, diameter_threshold=8, connectivity=1, parent=None, tree_traverser=None
  420. ):
  421. """Perform a diameter closing of the image.
  422. Diameter closing removes all dark structures of an image with
  423. maximal extension smaller than diameter_threshold. The maximal
  424. extension is defined as the maximal extension of the bounding box.
  425. The operator is also called Bounding Box Closing. In practice,
  426. the result is similar to a morphological closing, but long and thin
  427. structures are not removed.
  428. Technically, this operator is based on the max-tree representation of
  429. the image.
  430. Parameters
  431. ----------
  432. image : ndarray
  433. The input image for which the diameter_closing is to be calculated.
  434. This image can be of any type.
  435. diameter_threshold : unsigned int
  436. The maximal extension parameter (number of pixels). The default value
  437. is 8.
  438. connectivity : unsigned int, optional
  439. The neighborhood connectivity. The integer represents the maximum
  440. number of orthogonal steps to reach a neighbor. In 2D, it is 1 for
  441. a 4-neighborhood and 2 for a 8-neighborhood. Default value is 1.
  442. parent : ndarray, int64, optional
  443. Precomputed parent image representing the max tree of the inverted
  444. image. This function is fast, if precomputed parent and tree_traverser
  445. are provided. See Note for further details.
  446. tree_traverser : 1D array, int64, optional
  447. Precomputed traverser, where the pixels are ordered such that every
  448. pixel is preceded by its parent (except for the root which has no
  449. parent). This function is fast, if precomputed parent and
  450. tree_traverser are provided. See Note for further details.
  451. Returns
  452. -------
  453. output : ndarray
  454. Output image of the same shape and type as input image.
  455. See Also
  456. --------
  457. skimage.morphology.area_opening
  458. skimage.morphology.area_closing
  459. skimage.morphology.diameter_opening
  460. skimage.morphology.max_tree
  461. References
  462. ----------
  463. .. [1] Walter, T., & Klein, J.-C. (2002). Automatic Detection of
  464. Microaneurysms in Color Fundus Images of the Human Retina by Means
  465. of the Bounding Box Closing. In A. Colosimo, P. Sirabella,
  466. A. Giuliani (Eds.), Medical Data Analysis. Lecture Notes in Computer
  467. Science, vol 2526, pp. 210-220. Springer Berlin Heidelberg.
  468. :DOI:`10.1007/3-540-36104-9_23`
  469. .. [2] Carlinet, E., & Geraud, T. (2014). A Comparative Review of
  470. Component Tree Computation Algorithms. IEEE Transactions on Image
  471. Processing, 23(9), 3885-3895.
  472. :DOI:`10.1109/TIP.2014.2336551`
  473. Examples
  474. --------
  475. We create an image (quadratic function with a minimum in the center and
  476. 4 additional local minima.
  477. >>> w = 12
  478. >>> x, y = np.mgrid[0:w,0:w]
  479. >>> f = 180 + 0.2*((x - w/2)**2 + (y-w/2)**2)
  480. >>> f[2:3,1:5] = 160; f[2:4,9:11] = 140; f[9:11,2:4] = 120
  481. >>> f[9:10,9:11] = 100; f[10,10] = 100
  482. >>> f = f.astype(int)
  483. We can calculate the diameter closing:
  484. >>> closed = diameter_closing(f, 3, connectivity=1)
  485. All small minima with a maximal extension of 2 or less are removed.
  486. The remaining minima have all a maximal extension of at least 3.
  487. Notes
  488. -----
  489. If a max-tree representation (parent and tree_traverser) are given to the
  490. function, they must be calculated from the inverted image for this
  491. function, i.e.:
  492. >>> P, S = max_tree(invert(f))
  493. >>> closed = diameter_closing(f, 3, parent=P, tree_traverser=S)
  494. """
  495. # inversion of the input image
  496. image_inv = invert(image)
  497. output = image_inv.copy()
  498. if parent is None or tree_traverser is None:
  499. parent, tree_traverser = max_tree(image_inv, connectivity)
  500. diam = _max_tree._compute_extension(
  501. image_inv.ravel(),
  502. np.array(image_inv.shape, dtype=np.int32),
  503. parent.ravel(),
  504. tree_traverser,
  505. )
  506. _max_tree._direct_filter(
  507. image_inv.ravel(),
  508. output.ravel(),
  509. parent.ravel(),
  510. tree_traverser,
  511. diam,
  512. diameter_threshold,
  513. )
  514. output = invert(output)
  515. return output
  516. def max_tree_local_maxima(image, connectivity=1, parent=None, tree_traverser=None):
  517. """Determine all local maxima of the image.
  518. The local maxima are defined as connected sets of pixels with equal
  519. gray level strictly greater than the gray levels of all pixels in direct
  520. neighborhood of the set. The function labels the local maxima.
  521. Technically, the implementation is based on the max-tree representation
  522. of an image. The function is very efficient if the max-tree representation
  523. has already been computed. Otherwise, it is preferable to use
  524. the function local_maxima.
  525. Parameters
  526. ----------
  527. image : ndarray
  528. The input image for which the maxima are to be calculated.
  529. connectivity : unsigned int, optional
  530. The neighborhood connectivity. The integer represents the maximum
  531. number of orthogonal steps to reach a neighbor. In 2D, it is 1 for
  532. a 4-neighborhood and 2 for a 8-neighborhood. Default value is 1.
  533. parent : ndarray, int64, optional
  534. The value of each pixel is the index of its parent in the ravelled
  535. array.
  536. tree_traverser : 1D array, int64, optional
  537. The ordered pixel indices (referring to the ravelled array). The pixels
  538. are ordered such that every pixel is preceded by its parent (except for
  539. the root which has no parent).
  540. Returns
  541. -------
  542. local_max : ndarray, uint64
  543. Labeled local maxima of the image.
  544. See Also
  545. --------
  546. skimage.morphology.local_maxima
  547. skimage.morphology.max_tree
  548. References
  549. ----------
  550. .. [1] Vincent L., Proc. "Grayscale area openings and closings,
  551. their efficient implementation and applications",
  552. EURASIP Workshop on Mathematical Morphology and its
  553. Applications to Signal Processing, Barcelona, Spain, pp.22-27,
  554. May 1993.
  555. .. [2] Soille, P., "Morphological Image Analysis: Principles and
  556. Applications" (Chapter 6), 2nd edition (2003), ISBN 3540429883.
  557. :DOI:`10.1007/978-3-662-05088-0`
  558. .. [3] Salembier, P., Oliveras, A., & Garrido, L. (1998). Antiextensive
  559. Connected Operators for Image and Sequence Processing.
  560. IEEE Transactions on Image Processing, 7(4), 555-570.
  561. :DOI:`10.1109/83.663500`
  562. .. [4] Najman, L., & Couprie, M. (2006). Building the component tree in
  563. quasi-linear time. IEEE Transactions on Image Processing, 15(11),
  564. 3531-3539.
  565. :DOI:`10.1109/TIP.2006.877518`
  566. .. [5] Carlinet, E., & Geraud, T. (2014). A Comparative Review of
  567. Component Tree Computation Algorithms. IEEE Transactions on Image
  568. Processing, 23(9), 3885-3895.
  569. :DOI:`10.1109/TIP.2014.2336551`
  570. Examples
  571. --------
  572. We create an image (quadratic function with a maximum in the center and
  573. 4 additional constant maxima.
  574. >>> w = 10
  575. >>> x, y = np.mgrid[0:w,0:w]
  576. >>> f = 20 - 0.2*((x - w/2)**2 + (y-w/2)**2)
  577. >>> f[2:4,2:4] = 40; f[2:4,7:9] = 60; f[7:9,2:4] = 80; f[7:9,7:9] = 100
  578. >>> f = f.astype(int)
  579. We can calculate all local maxima:
  580. >>> maxima = max_tree_local_maxima(f)
  581. The resulting image contains the labeled local maxima.
  582. """
  583. output = np.ones(image.shape, dtype=np.uint64)
  584. if parent is None or tree_traverser is None:
  585. parent, tree_traverser = max_tree(image, connectivity)
  586. _max_tree._max_tree_local_maxima(
  587. image.ravel(), output.ravel(), parent.ravel(), tree_traverser
  588. )
  589. return output