_canny.py 9.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262
  1. """
  2. canny.py - Canny Edge detector
  3. Reference: Canny, J., A Computational Approach To Edge Detection, IEEE Trans.
  4. Pattern Analysis and Machine Intelligence, 8:679-714, 1986
  5. """
  6. import numpy as np
  7. import scipy.ndimage as ndi
  8. from ..util.dtype import dtype_limits
  9. from .._shared.filters import gaussian
  10. from .._shared.utils import _supported_float_type, check_nD
  11. from ._canny_cy import _nonmaximum_suppression_bilinear
  12. def _preprocess(image, mask, sigma, mode, cval):
  13. """Generate a smoothed image and an eroded mask.
  14. The image is smoothed using a gaussian filter ignoring masked
  15. pixels and the mask is eroded.
  16. Parameters
  17. ----------
  18. image : array
  19. Image to be smoothed.
  20. mask : array
  21. Mask with 1's for significant pixels, 0's for masked pixels.
  22. sigma : scalar or sequence of scalars
  23. Standard deviation for Gaussian kernel. The standard
  24. deviations of the Gaussian filter are given for each axis as a
  25. sequence, or as a single number, in which case it is equal for
  26. all axes.
  27. mode : str, {'reflect', 'constant', 'nearest', 'mirror', 'wrap'}
  28. The ``mode`` parameter determines how the array borders are
  29. handled, where ``cval`` is the value when mode is equal to
  30. 'constant'.
  31. cval : float, optional
  32. Value to fill past edges of input if `mode` is 'constant'.
  33. Returns
  34. -------
  35. smoothed_image : ndarray
  36. The smoothed array
  37. eroded_mask : ndarray
  38. The eroded mask.
  39. Notes
  40. -----
  41. This function calculates the fractional contribution of masked pixels
  42. by applying the function to the mask (which gets you the fraction of
  43. the pixel data that's due to significant points). We then mask the image
  44. and apply the function. The resulting values will be lower by the
  45. bleed-over fraction, so you can recalibrate by dividing by the function
  46. on the mask to recover the effect of smoothing from just the significant
  47. pixels.
  48. """
  49. gaussian_kwargs = dict(sigma=sigma, mode=mode, cval=cval, preserve_range=False)
  50. compute_bleedover = mode == 'constant' or mask is not None
  51. float_type = _supported_float_type(image.dtype)
  52. if mask is None:
  53. if compute_bleedover:
  54. mask = np.ones(image.shape, dtype=float_type)
  55. masked_image = image
  56. eroded_mask = np.ones(image.shape, dtype=bool)
  57. eroded_mask[:1, :] = 0
  58. eroded_mask[-1:, :] = 0
  59. eroded_mask[:, :1] = 0
  60. eroded_mask[:, -1:] = 0
  61. else:
  62. mask = mask.astype(bool, copy=False)
  63. masked_image = np.zeros_like(image)
  64. masked_image[mask] = image[mask]
  65. # Make the eroded mask. Setting the border value to zero will wipe
  66. # out the image edges for us.
  67. s = ndi.generate_binary_structure(2, 2)
  68. eroded_mask = ndi.binary_erosion(mask, s, border_value=0)
  69. if compute_bleedover:
  70. # Compute the fractional contribution of masked pixels by applying
  71. # the function to the mask (which gets you the fraction of the
  72. # pixel data that's due to significant points)
  73. bleed_over = (
  74. gaussian(mask.astype(float_type, copy=False), **gaussian_kwargs)
  75. + np.finfo(float_type).eps
  76. )
  77. # Smooth the masked image
  78. smoothed_image = gaussian(masked_image, **gaussian_kwargs)
  79. # Lower the result by the bleed-over fraction, so you can
  80. # recalibrate by dividing by the function on the mask to recover
  81. # the effect of smoothing from just the significant pixels.
  82. if compute_bleedover:
  83. smoothed_image /= bleed_over
  84. return smoothed_image, eroded_mask
  85. def canny(
  86. image,
  87. sigma=1.0,
  88. low_threshold=None,
  89. high_threshold=None,
  90. mask=None,
  91. use_quantiles=False,
  92. *,
  93. mode='constant',
  94. cval=0.0,
  95. ):
  96. """Edge filter an image using the Canny algorithm.
  97. Parameters
  98. ----------
  99. image : 2D array
  100. Grayscale input image to detect edges on; can be of any dtype.
  101. sigma : float, optional
  102. Standard deviation of the Gaussian filter.
  103. low_threshold : float, optional
  104. Lower bound for hysteresis thresholding (linking edges).
  105. If None, low_threshold is set to 10% of dtype's max.
  106. high_threshold : float, optional
  107. Upper bound for hysteresis thresholding (linking edges).
  108. If None, high_threshold is set to 20% of dtype's max.
  109. mask : array, dtype=bool, optional
  110. Mask to limit the application of Canny to a certain area.
  111. use_quantiles : bool, optional
  112. If ``True`` then treat low_threshold and high_threshold as
  113. quantiles of the edge magnitude image, rather than absolute
  114. edge magnitude values. If ``True`` then the thresholds must be
  115. in the range [0, 1].
  116. mode : str, {'reflect', 'constant', 'nearest', 'mirror', 'wrap'}
  117. The ``mode`` parameter determines how the array borders are
  118. handled during Gaussian filtering, where ``cval`` is the value when
  119. mode is equal to 'constant'.
  120. cval : float, optional
  121. Value to fill past edges of input if `mode` is 'constant'.
  122. Returns
  123. -------
  124. output : 2D array (image)
  125. The binary edge map.
  126. See also
  127. --------
  128. skimage.filters.sobel
  129. Notes
  130. -----
  131. The steps of the algorithm are as follows:
  132. * Smooth the image using a Gaussian with ``sigma`` width.
  133. * Apply the horizontal and vertical Sobel operators to get the gradients
  134. within the image. The edge strength is the norm of the gradient.
  135. * Thin potential edges to 1-pixel wide curves. First, find the normal
  136. to the edge at each point. This is done by looking at the
  137. signs and the relative magnitude of the X-Sobel and Y-Sobel
  138. to sort the points into 4 categories: horizontal, vertical,
  139. diagonal and antidiagonal. Then look in the normal and reverse
  140. directions to see if the values in either of those directions are
  141. greater than the point in question. Use interpolation to get a mix of
  142. points instead of picking the one that's the closest to the normal.
  143. * Perform a hysteresis thresholding: first label all points above the
  144. high threshold as edges. Then recursively label any point above the
  145. low threshold that is 8-connected to a labeled point as an edge.
  146. References
  147. ----------
  148. .. [1] Canny, J., A Computational Approach To Edge Detection, IEEE Trans.
  149. Pattern Analysis and Machine Intelligence, 8:679-714, 1986
  150. :DOI:`10.1109/TPAMI.1986.4767851`
  151. .. [2] William Green's Canny tutorial
  152. https://en.wikipedia.org/wiki/Canny_edge_detector
  153. Examples
  154. --------
  155. >>> from skimage import feature
  156. >>> rng = np.random.default_rng()
  157. >>> # Generate noisy image of a square
  158. >>> im = np.zeros((256, 256))
  159. >>> im[64:-64, 64:-64] = 1
  160. >>> im += 0.2 * rng.random(im.shape)
  161. >>> # First trial with the Canny filter, with the default smoothing
  162. >>> edges1 = feature.canny(im)
  163. >>> # Increase the smoothing for better results
  164. >>> edges2 = feature.canny(im, sigma=3)
  165. """
  166. # Regarding masks, any point touching a masked point will have a gradient
  167. # that is "infected" by the masked point, so it's enough to erode the
  168. # mask by one and then mask the output. We also mask out the border points
  169. # because who knows what lies beyond the edge of the image?
  170. if np.issubdtype(image.dtype, np.int64) or np.issubdtype(image.dtype, np.uint64):
  171. raise ValueError("64-bit integer images are not supported")
  172. check_nD(image, 2)
  173. dtype_max = dtype_limits(image, clip_negative=False)[1]
  174. if low_threshold is None:
  175. low_threshold = 0.1
  176. elif use_quantiles:
  177. if not (0.0 <= low_threshold <= 1.0):
  178. raise ValueError("Quantile thresholds must be between 0 and 1.")
  179. else:
  180. low_threshold /= dtype_max
  181. if high_threshold is None:
  182. high_threshold = 0.2
  183. elif use_quantiles:
  184. if not (0.0 <= high_threshold <= 1.0):
  185. raise ValueError("Quantile thresholds must be between 0 and 1.")
  186. else:
  187. high_threshold /= dtype_max
  188. if high_threshold < low_threshold:
  189. raise ValueError("low_threshold should be lower then high_threshold")
  190. # Image filtering
  191. smoothed, eroded_mask = _preprocess(image, mask, sigma, mode, cval)
  192. # Gradient magnitude estimation
  193. jsobel = ndi.sobel(smoothed, axis=1)
  194. isobel = ndi.sobel(smoothed, axis=0)
  195. magnitude = isobel * isobel
  196. magnitude += jsobel * jsobel
  197. np.sqrt(magnitude, out=magnitude)
  198. if use_quantiles:
  199. low_threshold, high_threshold = np.percentile(
  200. magnitude, [100.0 * low_threshold, 100.0 * high_threshold]
  201. )
  202. # Non-maximum suppression
  203. low_masked = _nonmaximum_suppression_bilinear(
  204. isobel, jsobel, magnitude, eroded_mask, low_threshold
  205. )
  206. # Double thresholding and edge tracking
  207. #
  208. # Segment the low-mask, then only keep low-segments that have
  209. # some high_mask component in them
  210. #
  211. low_mask = low_masked > 0
  212. strel = np.ones((3, 3), bool)
  213. labels, count = ndi.label(low_mask, strel)
  214. if count == 0:
  215. return low_mask
  216. high_mask = low_mask & (low_masked >= high_threshold)
  217. nonzero_sums = np.unique(labels[high_mask])
  218. good_label = np.zeros((count + 1,), bool)
  219. good_label[nonzero_sums] = True
  220. output_mask = good_label[labels]
  221. return output_mask