_util.py 12 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328
  1. """Utility functions used in the morphology subpackage."""
  2. import numpy as np
  3. from scipy import ndimage as ndi
  4. def _validate_connectivity(image_dim, connectivity, offset):
  5. """Convert any valid connectivity to a footprint and offset.
  6. Parameters
  7. ----------
  8. image_dim : int
  9. The number of dimensions of the input image.
  10. connectivity : int, array, or None
  11. The neighborhood connectivity. An integer is interpreted as in
  12. ``scipy.ndimage.generate_binary_structure``, as the maximum number
  13. of orthogonal steps to reach a neighbor. An array is directly
  14. interpreted as a footprint and its shape is validated against
  15. the input image shape. ``None`` is interpreted as a connectivity of 1.
  16. offset : tuple of int, or None
  17. The coordinates of the center of the footprint.
  18. Returns
  19. -------
  20. c_connectivity : array of bool
  21. The footprint (structuring element) corresponding to the input
  22. `connectivity`.
  23. offset : array of int
  24. The offset corresponding to the center of the footprint.
  25. Raises
  26. ------
  27. ValueError:
  28. If the image dimension and the connectivity or offset dimensions don't
  29. match.
  30. """
  31. if connectivity is None:
  32. connectivity = 1
  33. if np.isscalar(connectivity):
  34. c_connectivity = ndi.generate_binary_structure(image_dim, connectivity)
  35. else:
  36. c_connectivity = np.array(connectivity, bool)
  37. if c_connectivity.ndim != image_dim:
  38. raise ValueError("Connectivity dimension must be same as image")
  39. if offset is None:
  40. if any([x % 2 == 0 for x in c_connectivity.shape]):
  41. raise ValueError("Connectivity array must have an unambiguous " "center")
  42. offset = np.array(c_connectivity.shape) // 2
  43. return c_connectivity, offset
  44. def _raveled_offsets_and_distances(
  45. image_shape,
  46. *,
  47. footprint=None,
  48. connectivity=1,
  49. center=None,
  50. spacing=None,
  51. order='C',
  52. ):
  53. """Compute offsets to neighboring pixels in raveled coordinate space.
  54. This function also returns the corresponding distances from the center
  55. pixel given a spacing (assumed to be 1 along each axis by default).
  56. Parameters
  57. ----------
  58. image_shape : tuple of int
  59. The shape of the image for which the offsets are being computed.
  60. footprint : array of bool
  61. The footprint of the neighborhood, expressed as an n-dimensional array
  62. of 1s and 0s. If provided, the connectivity argument is ignored.
  63. connectivity : {1, ..., ndim}
  64. The square connectivity of the neighborhood: the number of orthogonal
  65. steps allowed to consider a pixel a neighbor. See
  66. `scipy.ndimage.generate_binary_structure`. Ignored if footprint is
  67. provided.
  68. center : tuple of int
  69. Tuple of indices to the center of the footprint. If not provided, it
  70. is assumed to be the center of the footprint, either provided or
  71. generated by the connectivity argument.
  72. spacing : tuple of float
  73. The spacing between pixels/voxels along each axis.
  74. order : 'C' or 'F'
  75. The ordering of the array, either C or Fortran ordering.
  76. Returns
  77. -------
  78. raveled_offsets : ndarray
  79. Linear offsets to a samples neighbors in the raveled image, sorted by
  80. their distance from the center.
  81. distances : ndarray
  82. The pixel distances corresponding to each offset.
  83. Notes
  84. -----
  85. This function will return values even if `image_shape` contains a dimension
  86. length that is smaller than `footprint`.
  87. Examples
  88. --------
  89. >>> off, d = _raveled_offsets_and_distances(
  90. ... (4, 5), footprint=np.ones((4, 3)), center=(1, 1)
  91. ... )
  92. >>> off
  93. array([-5, -1, 1, 5, -6, -4, 4, 6, 10, 9, 11])
  94. >>> d[0]
  95. 1.0
  96. >>> d[-1] # distance from (1, 1) to (3, 2)
  97. 2.236...
  98. """
  99. ndim = len(image_shape)
  100. if footprint is None:
  101. footprint = ndi.generate_binary_structure(rank=ndim, connectivity=connectivity)
  102. if center is None:
  103. center = tuple(s // 2 for s in footprint.shape)
  104. if not footprint.ndim == ndim == len(center):
  105. raise ValueError(
  106. "number of dimensions in image shape, footprint and its"
  107. "center index does not match"
  108. )
  109. offsets = np.stack(
  110. [(idx - c) for idx, c in zip(np.nonzero(footprint), center)], axis=-1
  111. )
  112. if order == 'F':
  113. offsets = offsets[:, ::-1]
  114. image_shape = image_shape[::-1]
  115. elif order != 'C':
  116. raise ValueError("order must be 'C' or 'F'")
  117. # Scale offsets in each dimension and sum
  118. ravel_factors = image_shape[1:] + (1,)
  119. ravel_factors = np.cumprod(ravel_factors[::-1])[::-1]
  120. raveled_offsets = (offsets * ravel_factors).sum(axis=1)
  121. # Sort by distance
  122. if spacing is None:
  123. spacing = np.ones(ndim)
  124. weighted_offsets = offsets * spacing
  125. distances = np.sqrt(np.sum(weighted_offsets**2, axis=1))
  126. sorted_raveled_offsets = raveled_offsets[np.argsort(distances, kind="stable")]
  127. sorted_distances = np.sort(distances, kind="stable")
  128. # If any dimension in image_shape is smaller than footprint.shape
  129. # duplicates might occur, remove them
  130. if any(x < y for x, y in zip(image_shape, footprint.shape)):
  131. # np.unique reorders, which we don't want
  132. _, indices = np.unique(sorted_raveled_offsets, return_index=True)
  133. indices = np.sort(indices, kind="stable")
  134. sorted_raveled_offsets = sorted_raveled_offsets[indices]
  135. sorted_distances = sorted_distances[indices]
  136. # Remove "offset to center"
  137. sorted_raveled_offsets = sorted_raveled_offsets[1:]
  138. sorted_distances = sorted_distances[1:]
  139. return sorted_raveled_offsets, sorted_distances
  140. def _offsets_to_raveled_neighbors(image_shape, footprint, center, order='C'):
  141. """Compute offsets to a samples neighbors if the image would be raveled.
  142. Parameters
  143. ----------
  144. image_shape : tuple
  145. The shape of the image for which the offsets are computed.
  146. footprint : ndarray
  147. The footprint (structuring element) determining the neighborhood
  148. expressed as an n-D array of 1's and 0's.
  149. center : tuple
  150. Tuple of indices to the center of `footprint`.
  151. order : {"C", "F"}, optional
  152. Whether the image described by `image_shape` is in row-major (C-style)
  153. or column-major (Fortran-style) order.
  154. Returns
  155. -------
  156. raveled_offsets : ndarray
  157. Linear offsets to a samples neighbors in the raveled image, sorted by
  158. their distance from the center.
  159. Notes
  160. -----
  161. This function will return values even if `image_shape` contains a dimension
  162. length that is smaller than `footprint`.
  163. Examples
  164. --------
  165. >>> _offsets_to_raveled_neighbors((4, 5), np.ones((4, 3)), (1, 1))
  166. array([-5, -1, 1, 5, -6, -4, 4, 6, 10, 9, 11])
  167. >>> _offsets_to_raveled_neighbors((2, 3, 2), np.ones((3, 3, 3)), (1, 1, 1))
  168. array([-6, -2, -1, 1, 2, 6, -8, -7, -5, -4, -3, 3, 4, 5, 7, 8, -9,
  169. 9])
  170. """
  171. raveled_offsets = _raveled_offsets_and_distances(
  172. image_shape, footprint=footprint, center=center, order=order
  173. )[0]
  174. return raveled_offsets
  175. def _resolve_neighborhood(footprint, connectivity, ndim, enforce_adjacency=True):
  176. """Validate or create a footprint (structuring element).
  177. Depending on the values of `connectivity` and `footprint` this function
  178. either creates a new footprint (`footprint` is None) using `connectivity`
  179. or validates the given footprint (`footprint` is not None).
  180. Parameters
  181. ----------
  182. footprint : ndarray
  183. The footprint (structuring) element used to determine the neighborhood
  184. of each evaluated pixel (``True`` denotes a connected pixel). It must
  185. be a boolean array and have the same number of dimensions as `image`.
  186. If neither `footprint` nor `connectivity` are given, all adjacent
  187. pixels are considered as part of the neighborhood.
  188. connectivity : int
  189. A number used to determine the neighborhood of each evaluated pixel.
  190. Adjacent pixels whose squared distance from the center is less than or
  191. equal to `connectivity` are considered neighbors. Ignored if
  192. `footprint` is not None.
  193. ndim : int
  194. Number of dimensions `footprint` ought to have.
  195. enforce_adjacency : bool
  196. A boolean that determines whether footprint must only specify direct
  197. neighbors.
  198. Returns
  199. -------
  200. footprint : ndarray
  201. Validated or new footprint specifying the neighborhood.
  202. Examples
  203. --------
  204. >>> _resolve_neighborhood(None, 1, 2)
  205. array([[False, True, False],
  206. [ True, True, True],
  207. [False, True, False]])
  208. >>> _resolve_neighborhood(None, None, 3).shape
  209. (3, 3, 3)
  210. """
  211. if footprint is None:
  212. if connectivity is None:
  213. connectivity = ndim
  214. footprint = ndi.generate_binary_structure(ndim, connectivity)
  215. else:
  216. # Validate custom structured element
  217. footprint = np.asarray(footprint, dtype=bool)
  218. # Must specify neighbors for all dimensions
  219. if footprint.ndim != ndim:
  220. raise ValueError(
  221. "number of dimensions in image and footprint do not" "match"
  222. )
  223. # Must only specify direct neighbors
  224. if enforce_adjacency and any(s != 3 for s in footprint.shape):
  225. raise ValueError("dimension size in footprint is not 3")
  226. elif any((s % 2 != 1) for s in footprint.shape):
  227. raise ValueError("footprint size must be odd along all dimensions")
  228. return footprint
  229. def _set_border_values(image, value, border_width=1):
  230. """Set edge values along all axes to a constant value.
  231. Parameters
  232. ----------
  233. image : ndarray
  234. The array to modify inplace.
  235. value : scalar
  236. The value to use. Should be compatible with `image`'s dtype.
  237. border_width : int or sequence of tuples
  238. A sequence with one 2-tuple per axis where the first and second values
  239. are the width of the border at the start and end of the axis,
  240. respectively. If an int is provided, a uniform border width along all
  241. axes is used.
  242. Examples
  243. --------
  244. >>> image = np.zeros((4, 5), dtype=int)
  245. >>> _set_border_values(image, 1)
  246. >>> image
  247. array([[1, 1, 1, 1, 1],
  248. [1, 0, 0, 0, 1],
  249. [1, 0, 0, 0, 1],
  250. [1, 1, 1, 1, 1]])
  251. >>> image = np.zeros((8, 8), dtype=int)
  252. >>> _set_border_values(image, 1, border_width=((1, 1), (2, 3)))
  253. >>> image
  254. array([[1, 1, 1, 1, 1, 1, 1, 1],
  255. [1, 1, 0, 0, 0, 1, 1, 1],
  256. [1, 1, 0, 0, 0, 1, 1, 1],
  257. [1, 1, 0, 0, 0, 1, 1, 1],
  258. [1, 1, 0, 0, 0, 1, 1, 1],
  259. [1, 1, 0, 0, 0, 1, 1, 1],
  260. [1, 1, 0, 0, 0, 1, 1, 1],
  261. [1, 1, 1, 1, 1, 1, 1, 1]])
  262. """
  263. if np.isscalar(border_width):
  264. border_width = ((border_width, border_width),) * image.ndim
  265. elif len(border_width) != image.ndim:
  266. raise ValueError('length of `border_width` must match image.ndim')
  267. for axis, npad in enumerate(border_width):
  268. if len(npad) != 2:
  269. raise ValueError('each sequence in `border_width` must have ' 'length 2')
  270. w_start, w_end = npad
  271. if w_start == w_end == 0:
  272. continue
  273. elif w_start == w_end == 1:
  274. # Index first and last element in the current dimension
  275. sl = (slice(None),) * axis + ((0, -1),) + (...,)
  276. image[sl] = value
  277. continue
  278. if w_start > 0:
  279. # set first w_start entries along axis to value
  280. sl = (slice(None),) * axis + (slice(0, w_start),) + (...,)
  281. image[sl] = value
  282. if w_end > 0:
  283. # set last w_end entries along axis to value
  284. sl = (slice(None),) * axis + (slice(-w_end, None),) + (...,)
  285. image[sl] = value