_basic_features.py 6.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198
  1. from itertools import combinations_with_replacement
  2. import itertools
  3. import numpy as np
  4. from skimage import filters, feature
  5. from skimage.util.dtype import img_as_float32
  6. from .._shared._dependency_checks import is_wasm
  7. if not is_wasm:
  8. from concurrent.futures import ThreadPoolExecutor as PoolExecutor
  9. else:
  10. from contextlib import AbstractContextManager
  11. # Threading isn't supported on WASM, mock ThreadPoolExecutor as a fallback
  12. class PoolExecutor(AbstractContextManager):
  13. def __init__(self, *_, **__):
  14. pass
  15. def __exit__(self, exc_type, exc_val, exc_tb):
  16. pass
  17. def map(self, fn, iterables):
  18. return map(fn, iterables)
  19. def _texture_filter(gaussian_filtered):
  20. H_elems = [
  21. np.gradient(np.gradient(gaussian_filtered)[ax0], axis=ax1)
  22. for ax0, ax1 in combinations_with_replacement(range(gaussian_filtered.ndim), 2)
  23. ]
  24. eigvals = feature.hessian_matrix_eigvals(H_elems)
  25. return eigvals
  26. def _singlescale_basic_features_singlechannel(
  27. img, sigma, intensity=True, edges=True, texture=True
  28. ):
  29. results = ()
  30. gaussian_filtered = filters.gaussian(img, sigma=sigma, preserve_range=False)
  31. if intensity:
  32. results += (gaussian_filtered,)
  33. if edges:
  34. results += (filters.sobel(gaussian_filtered),)
  35. if texture:
  36. results += (*_texture_filter(gaussian_filtered),)
  37. return results
  38. def _mutiscale_basic_features_singlechannel(
  39. img,
  40. intensity=True,
  41. edges=True,
  42. texture=True,
  43. sigma_min=0.5,
  44. sigma_max=16,
  45. num_sigma=None,
  46. workers=None,
  47. ):
  48. """Features for a single channel nd image.
  49. Parameters
  50. ----------
  51. img : ndarray
  52. Input image, which can be grayscale or multichannel.
  53. intensity : bool, default True
  54. If True, pixel intensities averaged over the different scales
  55. are added to the feature set.
  56. edges : bool, default True
  57. If True, intensities of local gradients averaged over the different
  58. scales are added to the feature set.
  59. texture : bool, default True
  60. If True, eigenvalues of the Hessian matrix after Gaussian blurring
  61. at different scales are added to the feature set.
  62. sigma_min : float, optional
  63. Smallest value of the Gaussian kernel used to average local
  64. neighborhoods before extracting features.
  65. sigma_max : float, optional
  66. Largest value of the Gaussian kernel used to average local
  67. neighborhoods before extracting features.
  68. num_sigma : int, optional
  69. Number of values of the Gaussian kernel between sigma_min and sigma_max.
  70. If None, sigma_min multiplied by powers of 2 are used.
  71. workers : int or None, optional
  72. The number of parallel threads to use. If set to ``None``, the full
  73. set of available cores are used.
  74. Returns
  75. -------
  76. features : list
  77. List of features, each element of the list is an array of shape as img.
  78. """
  79. # computations are faster as float32
  80. img = np.ascontiguousarray(img_as_float32(img))
  81. if num_sigma is None:
  82. num_sigma = int(np.log2(sigma_max) - np.log2(sigma_min) + 1)
  83. sigmas = np.logspace(
  84. np.log2(sigma_min),
  85. np.log2(sigma_max),
  86. num=num_sigma,
  87. base=2,
  88. endpoint=True,
  89. )
  90. with PoolExecutor(max_workers=workers) as ex:
  91. out_sigmas = list(
  92. ex.map(
  93. lambda s: _singlescale_basic_features_singlechannel(
  94. img, s, intensity=intensity, edges=edges, texture=texture
  95. ),
  96. sigmas,
  97. )
  98. )
  99. features = itertools.chain.from_iterable(out_sigmas)
  100. return features
  101. def multiscale_basic_features(
  102. image,
  103. intensity=True,
  104. edges=True,
  105. texture=True,
  106. sigma_min=0.5,
  107. sigma_max=16,
  108. num_sigma=None,
  109. workers=None,
  110. *,
  111. channel_axis=None,
  112. ):
  113. """Local features for a single- or multi-channel nd image.
  114. Intensity, gradient intensity and local structure are computed at
  115. different scales thanks to Gaussian blurring.
  116. Parameters
  117. ----------
  118. image : ndarray
  119. Input image, which can be grayscale or multichannel.
  120. intensity : bool, default True
  121. If True, pixel intensities averaged over the different scales
  122. are added to the feature set.
  123. edges : bool, default True
  124. If True, intensities of local gradients averaged over the different
  125. scales are added to the feature set.
  126. texture : bool, default True
  127. If True, eigenvalues of the Hessian matrix after Gaussian blurring
  128. at different scales are added to the feature set.
  129. sigma_min : float, optional
  130. Smallest value of the Gaussian kernel used to average local
  131. neighborhoods before extracting features.
  132. sigma_max : float, optional
  133. Largest value of the Gaussian kernel used to average local
  134. neighborhoods before extracting features.
  135. num_sigma : int, optional
  136. Number of values of the Gaussian kernel between sigma_min and sigma_max.
  137. If None, sigma_min multiplied by powers of 2 are used.
  138. workers : int or None, optional
  139. The number of parallel threads to use. If set to ``None``, the full
  140. set of available cores are used.
  141. channel_axis : int or None, optional
  142. If None, the image is assumed to be a grayscale (single channel) image.
  143. Otherwise, this parameter indicates which axis of the array corresponds
  144. to channels.
  145. .. versionadded:: 0.19
  146. ``channel_axis`` was added in 0.19.
  147. Returns
  148. -------
  149. features : np.ndarray
  150. Array of shape ``image.shape + (n_features,)``. When `channel_axis` is
  151. not None, all channels are concatenated along the features dimension.
  152. (i.e. ``n_features == n_features_singlechannel * n_channels``)
  153. """
  154. if not any([intensity, edges, texture]):
  155. raise ValueError(
  156. "At least one of `intensity`, `edges` or `textures`"
  157. "must be True for features to be computed."
  158. )
  159. if channel_axis is None:
  160. image = image[..., np.newaxis]
  161. channel_axis = -1
  162. elif channel_axis != -1:
  163. image = np.moveaxis(image, channel_axis, -1)
  164. all_results = (
  165. _mutiscale_basic_features_singlechannel(
  166. image[..., dim],
  167. intensity=intensity,
  168. edges=edges,
  169. texture=texture,
  170. sigma_min=sigma_min,
  171. sigma_max=sigma_max,
  172. num_sigma=num_sigma,
  173. workers=workers,
  174. )
  175. for dim in range(image.shape[-1])
  176. )
  177. features = list(itertools.chain.from_iterable(all_results))
  178. out = np.stack(features, axis=-1)
  179. return out