convex_hull.py 8.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222
  1. """Convex Hull."""
  2. from itertools import product
  3. import numpy as np
  4. from scipy.spatial import ConvexHull, QhullError
  5. from ..measure.pnpoly import grid_points_in_poly
  6. from ._convex_hull import possible_hull
  7. from ..measure._label import label
  8. from ..util import unique_rows
  9. from .._shared.utils import warn
  10. __all__ = ['convex_hull_image', 'convex_hull_object']
  11. def _offsets_diamond(ndim):
  12. offsets = np.zeros((2 * ndim, ndim))
  13. for vertex, (axis, offset) in enumerate(product(range(ndim), (-0.5, 0.5))):
  14. offsets[vertex, axis] = offset
  15. return offsets
  16. def _check_coords_in_hull(gridcoords, hull_equations, tolerance):
  17. r"""Checks all the coordinates for inclusiveness in the convex hull.
  18. Parameters
  19. ----------
  20. gridcoords : (M, N) ndarray
  21. Coordinates of ``N`` points in ``M`` dimensions.
  22. hull_equations : (M, N) ndarray
  23. Hyperplane equations of the facets of the convex hull.
  24. tolerance : float
  25. Tolerance when determining whether a point is inside the hull. Due
  26. to numerical floating point errors, a tolerance of 0 can result in
  27. some points erroneously being classified as being outside the hull.
  28. Returns
  29. -------
  30. coords_in_hull : ndarray of bool
  31. Binary 1D ndarray representing points in n-dimensional space
  32. with value ``True`` set for points inside the convex hull.
  33. Notes
  34. -----
  35. Checking the inclusiveness of coordinates in a convex hull requires
  36. intermediate calculations of dot products which are memory-intensive.
  37. Thus, the convex hull equations are checked individually with all
  38. coordinates to keep within the memory limit.
  39. References
  40. ----------
  41. .. [1] https://github.com/scikit-image/scikit-image/issues/5019
  42. """
  43. ndim, n_coords = gridcoords.shape
  44. n_hull_equations = hull_equations.shape[0]
  45. coords_in_hull = np.ones(n_coords, dtype=bool)
  46. # Pre-allocate arrays to cache intermediate results for reducing overheads
  47. dot_array = np.empty(n_coords, dtype=np.float64)
  48. test_ineq_temp = np.empty(n_coords, dtype=np.float64)
  49. coords_single_ineq = np.empty(n_coords, dtype=bool)
  50. # A point is in the hull if it satisfies all of the hull's inequalities
  51. for idx in range(n_hull_equations):
  52. # Tests a hyperplane equation on all coordinates of volume
  53. np.dot(hull_equations[idx, :ndim], gridcoords, out=dot_array)
  54. np.add(dot_array, hull_equations[idx, ndim:], out=test_ineq_temp)
  55. np.less(test_ineq_temp, tolerance, out=coords_single_ineq)
  56. coords_in_hull *= coords_single_ineq
  57. return coords_in_hull
  58. def convex_hull_image(
  59. image, offset_coordinates=True, tolerance=1e-10, include_borders=True
  60. ):
  61. """Compute the convex hull image of a binary image.
  62. The convex hull is the set of pixels included in the smallest convex
  63. polygon that surround all white pixels in the input image.
  64. Parameters
  65. ----------
  66. image : array
  67. Binary input image. This array is cast to bool before processing.
  68. offset_coordinates : bool, optional
  69. If ``True``, a pixel at coordinate, e.g., (4, 7) will be represented
  70. by coordinates (3.5, 7), (4.5, 7), (4, 6.5), and (4, 7.5). This adds
  71. some "extent" to a pixel when computing the hull.
  72. tolerance : float, optional
  73. Tolerance when determining whether a point is inside the hull. Due
  74. to numerical floating point errors, a tolerance of 0 can result in
  75. some points erroneously being classified as being outside the hull.
  76. include_borders : bool, optional
  77. If ``False``, vertices/edges are excluded from the final hull mask.
  78. Returns
  79. -------
  80. hull : (M, N) array of bool
  81. Binary image with pixels in convex hull set to True.
  82. References
  83. ----------
  84. .. [1] https://blogs.mathworks.com/steve/2011/10/04/binary-image-convex-hull-algorithm-notes/
  85. """
  86. ndim = image.ndim
  87. if np.count_nonzero(image) == 0:
  88. warn(
  89. "Input image is entirely zero, no valid convex hull. "
  90. "Returning empty image",
  91. UserWarning,
  92. )
  93. return np.zeros(image.shape, dtype=bool)
  94. # In 2D, we do an optimisation by choosing only pixels that are
  95. # the starting or ending pixel of a row or column. This vastly
  96. # limits the number of coordinates to examine for the virtual hull.
  97. if ndim == 2:
  98. coords = possible_hull(np.ascontiguousarray(image, dtype=np.uint8))
  99. else:
  100. coords = np.transpose(np.nonzero(image))
  101. if offset_coordinates:
  102. # when offsetting, we multiply number of vertices by 2 * ndim.
  103. # therefore, we reduce the number of coordinates by using a
  104. # convex hull on the original set, before offsetting.
  105. try:
  106. hull0 = ConvexHull(coords)
  107. except QhullError as err:
  108. warn(
  109. f"Failed to get convex hull image. "
  110. f"Returning empty image, see error message below:\n"
  111. f"{err}"
  112. )
  113. return np.zeros(image.shape, dtype=bool)
  114. coords = hull0.points[hull0.vertices]
  115. # Add a vertex for the middle of each pixel edge
  116. if offset_coordinates:
  117. offsets = _offsets_diamond(image.ndim)
  118. coords = (coords[:, np.newaxis, :] + offsets).reshape(-1, ndim)
  119. # repeated coordinates can *sometimes* cause problems in
  120. # scipy.spatial.ConvexHull, so we remove them.
  121. coords = unique_rows(coords)
  122. # Find the convex hull
  123. try:
  124. hull = ConvexHull(coords)
  125. except QhullError as err:
  126. warn(
  127. f"Failed to get convex hull image. "
  128. f"Returning empty image, see error message below:\n"
  129. f"{err}"
  130. )
  131. return np.zeros(image.shape, dtype=bool)
  132. vertices = hull.points[hull.vertices]
  133. # If 2D, use fast Cython function to locate convex hull pixels
  134. if ndim == 2:
  135. labels = grid_points_in_poly(image.shape, vertices, binarize=False)
  136. # If include_borders is True, we include vertices (2) and edge
  137. # points (3) in the mask, otherwise only the inside of the hull (1)
  138. mask = labels >= 1 if include_borders else labels == 1
  139. else:
  140. gridcoords = np.reshape(np.mgrid[tuple(map(slice, image.shape))], (ndim, -1))
  141. coords_in_hull = _check_coords_in_hull(gridcoords, hull.equations, tolerance)
  142. mask = np.reshape(coords_in_hull, image.shape)
  143. return mask
  144. def convex_hull_object(image, *, connectivity=2):
  145. r"""Compute the convex hull image of individual objects in a binary image.
  146. The convex hull is the set of pixels included in the smallest convex
  147. polygon that surround all white pixels in the input image.
  148. Parameters
  149. ----------
  150. image : (M, N) ndarray
  151. Binary input image.
  152. connectivity : {1, 2}, int, optional
  153. Determines the neighbors of each pixel. Adjacent elements
  154. within a squared distance of ``connectivity`` from pixel center
  155. are considered neighbors.::
  156. 1-connectivity 2-connectivity
  157. [ ] [ ] [ ] [ ]
  158. | \ | /
  159. [ ]--[x]--[ ] [ ]--[x]--[ ]
  160. | / | \
  161. [ ] [ ] [ ] [ ]
  162. Returns
  163. -------
  164. hull : ndarray of bool
  165. Binary image with pixels inside convex hull set to ``True``.
  166. Notes
  167. -----
  168. This function uses ``skimage.morphology.label`` to define unique objects,
  169. finds the convex hull of each using ``convex_hull_image``, and combines
  170. these regions with logical OR. Be aware the convex hulls of unconnected
  171. objects may overlap in the result. If this is suspected, consider using
  172. convex_hull_image separately on each object or adjust ``connectivity``.
  173. """
  174. if image.ndim > 2:
  175. raise ValueError("Input must be a 2D image")
  176. if connectivity not in (1, 2):
  177. raise ValueError('`connectivity` must be either 1 or 2.')
  178. labeled_im = label(image, connectivity=connectivity, background=0)
  179. convex_obj = np.zeros(image.shape, dtype=bool)
  180. convex_img = np.zeros(image.shape, dtype=bool)
  181. for i in range(1, labeled_im.max() + 1):
  182. convex_obj = convex_hull_image(labeled_im == i)
  183. convex_img = np.logical_or(convex_img, convex_obj)
  184. return convex_img