_rolling_ball.py 6.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211
  1. import numpy as np
  2. from .._shared.utils import _supported_float_type, deprecate_parameter, DEPRECATED
  3. from ._rolling_ball_cy import apply_kernel, apply_kernel_nan
  4. @deprecate_parameter(
  5. "num_threads", new_name="workers", start_version="0.26", stop_version="0.28"
  6. )
  7. def rolling_ball(
  8. image,
  9. *,
  10. radius=100,
  11. kernel=None,
  12. nansafe=False,
  13. num_threads=DEPRECATED,
  14. workers=None,
  15. ):
  16. """Estimate background intensity using the rolling-ball algorithm.
  17. This function is a generalization of the rolling-ball algorithm [1]_ to
  18. estimate the background intensity of an n-dimensional image. This is
  19. typically useful for background subtraction in case of uneven exposure.
  20. Think of the image as a landscape (where altitude is determined by
  21. intensity), under which a ball of given radius is rolled. At each
  22. position, the ball's apex gives the resulting background intensity.
  23. Parameters
  24. ----------
  25. image : ndarray
  26. The image to be filtered.
  27. radius : int, optional
  28. Radius of the ball-shaped kernel to be rolled under the
  29. image landscape. Used only if `kernel` is ``None``.
  30. kernel : ndarray, optional
  31. An alternative way to specify the rolling ball, as an arbitrary
  32. kernel. It must have the same number of axes as `image`.
  33. nansafe : bool, optional
  34. If ``False`` (default), the function assumes that none of the values
  35. in `image` are ``np.nan``, and uses a faster implementation.
  36. workers : int, optional
  37. The maximum number of threads to use. If ``None``, use the OpenMP
  38. default value; typically equal to the maximum number of virtual cores.
  39. Note: This is an upper limit to the number of threads. The exact number
  40. is determined by the system's OpenMP library.
  41. .. versionadded:: 0.26
  42. Replaces deprecated parameter `num_threads`.
  43. Returns
  44. -------
  45. background : ndarray
  46. The estimated background of the image.
  47. Notes
  48. -----
  49. This implementation assumes that dark pixels correspond to the background. If
  50. you have a bright background, invert the image before passing it to this
  51. function, e.g., using :func:`skimage.util.invert`.
  52. For this method to give meaningful results, the radius of the ball (or
  53. typical size of the kernel, in the general case) should be larger than the
  54. typical size of the image features of interest.
  55. This algorithm is sensitive to noise (in particular salt-and-pepper
  56. noise). If this is a problem in your image, you can apply mild
  57. Gaussian smoothing before passing the image to this function.
  58. This algorithm's complexity is polynomial in the radius, with degree equal
  59. to the image dimensionality (a 2D image is N^2, a 3D image is N^3, etc.),
  60. so it can take a long time as the radius grows beyond 30 or so ([2]_, [3]_).
  61. It is an exact N-dimensional calculation; if all you need is an
  62. approximation, faster options to consider are top-hat filtering [4]_ or
  63. downscaling-then-upscaling to reduce the size of the input processed.
  64. References
  65. ----------
  66. .. [1] Sternberg, Stanley R. "Biomedical image processing." Computer 1
  67. (1983): 22-34. :DOI:`10.1109/MC.1983.1654163`
  68. .. [2] https://github.com/scikit-image/scikit-image/issues/5193
  69. .. [3] https://github.com/scikit-image/scikit-image/issues/7423
  70. .. [4] https://forum.image.sc/t/59267/7
  71. Examples
  72. --------
  73. >>> import numpy as np
  74. >>> import skimage as ski
  75. >>> image = ski.data.coins()
  76. >>> background = ski.restoration.rolling_ball(image)
  77. >>> filtered_image = image - background
  78. >>> import numpy as np
  79. >>> import skimage as ski
  80. >>> image = ski.data.coins()
  81. >>> kernel = ski.restoration.ellipsoid_kernel((101, 101), 75)
  82. >>> background = ski.restoration.rolling_ball(image, kernel=kernel)
  83. >>> filtered_image = image - background
  84. """
  85. image = np.asarray(image)
  86. float_type = _supported_float_type(image.dtype)
  87. img = image.astype(float_type, copy=False)
  88. if workers is None:
  89. workers = 0
  90. if kernel is None:
  91. kernel = ball_kernel(radius, image.ndim)
  92. kernel = kernel.astype(float_type)
  93. kernel_shape = np.asarray(kernel.shape)
  94. kernel_center = kernel_shape // 2
  95. center_intensity = kernel[tuple(kernel_center)]
  96. intensity_difference = center_intensity - kernel
  97. intensity_difference[kernel == np.inf] = np.inf
  98. intensity_difference = intensity_difference.astype(img.dtype)
  99. intensity_difference = intensity_difference.reshape(-1)
  100. img = np.pad(
  101. img, kernel_center[:, np.newaxis], constant_values=np.inf, mode="constant"
  102. )
  103. func = apply_kernel_nan if nansafe else apply_kernel
  104. background = func(
  105. img.reshape(-1),
  106. intensity_difference,
  107. np.zeros_like(image, dtype=img.dtype).reshape(-1),
  108. np.array(image.shape, dtype=np.intp),
  109. np.array(img.shape, dtype=np.intp),
  110. kernel_shape.astype(np.intp),
  111. workers,
  112. )
  113. background = background.astype(image.dtype, copy=False)
  114. return background
  115. def ball_kernel(radius, ndim):
  116. """Create a ball kernel for restoration.rolling_ball.
  117. Parameters
  118. ----------
  119. radius : int
  120. Radius of the ball.
  121. ndim : int
  122. Number of dimensions of the ball. ``ndim`` should match the
  123. dimensionality of the image the kernel will be applied to.
  124. Returns
  125. -------
  126. kernel : ndarray
  127. The kernel containing the surface intensity of the top half
  128. of the ellipsoid.
  129. See Also
  130. --------
  131. rolling_ball
  132. """
  133. kernel_coords = np.stack(
  134. np.meshgrid(
  135. *[np.arange(-x, x + 1) for x in [np.ceil(radius)] * ndim], indexing='ij'
  136. ),
  137. axis=-1,
  138. )
  139. sum_of_squares = np.sum(kernel_coords**2, axis=-1)
  140. distance_from_center = np.sqrt(sum_of_squares)
  141. kernel = np.sqrt(np.clip(radius**2 - sum_of_squares, 0, None))
  142. kernel[distance_from_center > radius] = np.inf
  143. return kernel
  144. def ellipsoid_kernel(shape, intensity):
  145. """Create an ellipoid kernel for restoration.rolling_ball.
  146. Parameters
  147. ----------
  148. shape : array-like
  149. Length of the principal axis of the ellipsoid (excluding
  150. the intensity axis). The kernel needs to have the same
  151. dimensionality as the image it will be applied to.
  152. intensity : int
  153. Length of the intensity axis of the ellipsoid.
  154. Returns
  155. -------
  156. kernel : ndarray
  157. The kernel containing the surface intensity of the top half
  158. of the ellipsoid.
  159. See Also
  160. --------
  161. rolling_ball
  162. """
  163. shape = np.asarray(shape)
  164. semi_axis = np.clip(shape // 2, 1, None)
  165. kernel_coords = np.stack(
  166. np.meshgrid(*[np.arange(-x, x + 1) for x in semi_axis], indexing='ij'), axis=-1
  167. )
  168. intensity_scaling = 1 - np.sum((kernel_coords / semi_axis) ** 2, axis=-1)
  169. kernel = intensity * np.sqrt(np.clip(intensity_scaling, 0, None))
  170. kernel[intensity_scaling < 0] = np.inf
  171. return kernel