"Fossies" - the Fresh Open Source Software Archive

Member "scikit-image-0.19.3/skimage/feature/_canny.py" (12 Jun 2022, 13114 Bytes) of package /linux/misc/scikit-image-0.19.3.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Python source code syntax highlighting (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file. For more information about "_canny.py" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 0.19.2_vs_0.19.3.

    1 """
    2 canny.py - Canny Edge detector
    3 
    4 Reference: Canny, J., A Computational Approach To Edge Detection, IEEE Trans.
    5     Pattern Analysis and Machine Intelligence, 8:679-714, 1986
    6 
    7 Originally part of CellProfiler, code licensed under both GPL and BSD licenses.
    8 Website: http://www.cellprofiler.org
    9 Copyright (c) 2003-2009 Massachusetts Institute of Technology
   10 Copyright (c) 2009-2011 Broad Institute
   11 All rights reserved.
   12 Original author: Lee Kamentsky
   13 """
   14 
   15 import numpy as np
   16 import scipy.ndimage as ndi
   17 
   18 from ..util.dtype import dtype_limits
   19 from .._shared.filters import gaussian
   20 from .._shared.utils import check_nD
   21 from ..filters._gaussian import gaussian
   22 
   23 
   24 def _preprocess(image, mask, sigma, mode, cval):
   25     """Generate a smoothed image and an eroded mask.
   26 
   27     The image is smoothed using a gaussian filter ignoring masked
   28     pixels and the mask is eroded.
   29 
   30     Parameters
   31     ----------
   32     image : array
   33         Image to be smoothed.
   34     mask : array
   35         Mask with 1's for significant pixels, 0's for masked pixels.
   36     sigma : scalar or sequence of scalars
   37         Standard deviation for Gaussian kernel. The standard
   38         deviations of the Gaussian filter are given for each axis as a
   39         sequence, or as a single number, in which case it is equal for
   40         all axes.
   41     mode : str, {'reflect', 'constant', 'nearest', 'mirror', 'wrap'}
   42         The ``mode`` parameter determines how the array borders are
   43         handled, where ``cval`` is the value when mode is equal to
   44         'constant'.
   45     cval : float, optional
   46         Value to fill past edges of input if `mode` is 'constant'.
   47 
   48     Returns
   49     -------
   50     smoothed_image : ndarray
   51         The smoothed array
   52     eroded_mask : ndarray
   53         The eroded mask.
   54 
   55     Notes
   56     -----
   57     This function calculates the fractional contribution of masked pixels
   58     by applying the function to the mask (which gets you the fraction of
   59     the pixel data that's due to significant points). We then mask the image
   60     and apply the function. The resulting values will be lower by the
   61     bleed-over fraction, so you can recalibrate by dividing by the function
   62     on the mask to recover the effect of smoothing from just the significant
   63     pixels.
   64     """
   65 
   66     gaussian_kwargs = dict(sigma=sigma, mode=mode, cval=cval,
   67                            preserve_range=False)
   68     if mask is None:
   69         mask = np.ones(image.shape)
   70         masked_image = image.copy()
   71 
   72         eroded_mask = np.ones(image.shape, dtype=bool)
   73         eroded_mask[:1, :] = 0
   74         eroded_mask[-1:, :] = 0
   75         eroded_mask[:, :1] = 0
   76         eroded_mask[:, -1:] = 0
   77 
   78     else:
   79         mask = mask.astype(bool, copy=False)
   80         masked_image = np.zeros_like(image)
   81         masked_image[mask] = image[mask]
   82 
   83         # Make the eroded mask. Setting the border value to zero will wipe
   84         # out the image edges for us.
   85         s = ndi.generate_binary_structure(2, 2)
   86         eroded_mask = ndi.binary_erosion(mask, s, border_value=0)
   87 
   88     # Compute the fractional contribution of masked pixels by applying
   89     # the function to the mask (which gets you the fraction of the
   90     # pixel data that's due to significant points)
   91     bleed_over = gaussian(mask.astype(float, copy=False),
   92                           **gaussian_kwargs) + np.finfo(float).eps
   93 
   94     # Smooth the masked image
   95     smoothed_image = gaussian(masked_image, **gaussian_kwargs)
   96 
   97     # Lower the result by the bleed-over fraction, so you can
   98     # recalibrate by dividing by the function on the mask to recover
   99     # the effect of smoothing from just the significant pixels.
  100     smoothed_image /= bleed_over
  101 
  102     return smoothed_image, eroded_mask
  103 
  104 
  105 def _set_local_maxima(magnitude, pts, w_num, w_denum, row_slices,
  106                       col_slices, out):
  107     """Get the magnitudes shifted left to make a matrix of the points to
  108     the right of pts. Similarly, shift left and down to get the points
  109     to the top right of pts.
  110     """
  111     r_0, r_1, r_2, r_3 = row_slices
  112     c_0, c_1, c_2, c_3 = col_slices
  113     c1 = magnitude[r_0, c_0][pts[r_1, c_1]]
  114     c2 = magnitude[r_2, c_2][pts[r_3, c_3]]
  115     m = magnitude[pts]
  116     w = w_num[pts] / w_denum[pts]
  117     c_plus = c2 * w + c1 * (1 - w) <= m
  118     c1 = magnitude[r_1, c_1][pts[r_0, c_0]]
  119     c2 = magnitude[r_3, c_3][pts[r_2, c_2]]
  120     c_minus = c2 * w + c1 * (1 - w) <= m
  121     out[pts] = c_plus & c_minus
  122 
  123     return out
  124 
  125 
  126 def _get_local_maxima(isobel, jsobel, magnitude, eroded_mask):
  127     """Edge thinning by non-maximum suppression.
  128 
  129     Finds the normal to the edge at each point using the arctangent of the
  130     ratio of the Y sobel over the X sobel - pragmatically, we can
  131     look at the signs of X and Y and the relative magnitude of X vs Y
  132     to sort the points into 4 categories: horizontal, vertical,
  133     diagonal and antidiagonal.
  134 
  135     Look in the normal and reverse directions to see if the values
  136     in either of those directions are greater than the point in question.
  137     Use interpolation (via _set_local_maxima) to get a mix of points
  138     instead of picking the one that's the closest to the normal.
  139     """
  140     abs_isobel = np.abs(isobel)
  141     abs_jsobel = np.abs(jsobel)
  142 
  143     eroded_mask = eroded_mask & (magnitude > 0)
  144 
  145     # Normals' orientations
  146     is_horizontal = eroded_mask & (abs_isobel >= abs_jsobel)
  147     is_vertical = eroded_mask & (abs_isobel <= abs_jsobel)
  148     is_up = (isobel >= 0)
  149     is_down = (isobel <= 0)
  150     is_right = (jsobel >= 0)
  151     is_left = (jsobel <= 0)
  152     #
  153     # --------- Find local maxima --------------
  154     #
  155     # Assign each point to have a normal of 0-45 degrees, 45-90 degrees,
  156     # 90-135 degrees and 135-180 degrees.
  157     #
  158     local_maxima = np.zeros(magnitude.shape, bool)
  159     # ----- 0 to 45 degrees ------
  160     # Mix diagonal and horizontal
  161     pts_plus = is_up & is_right
  162     pts_minus = is_down & is_left
  163     pts = ((pts_plus | pts_minus) & is_horizontal)
  164     # Get the magnitudes shifted left to make a matrix of the points to the
  165     # right of pts. Similarly, shift left and down to get the points to the
  166     # top right of pts.
  167     local_maxima = _set_local_maxima(
  168         magnitude, pts, abs_jsobel, abs_isobel,
  169         [slice(1, None), slice(-1), slice(1, None), slice(-1)],
  170         [slice(None), slice(None), slice(1, None), slice(-1)],
  171         local_maxima)
  172     # ----- 45 to 90 degrees ------
  173     # Mix diagonal and vertical
  174     #
  175     pts = ((pts_plus | pts_minus) & is_vertical)
  176     local_maxima = _set_local_maxima(
  177         magnitude, pts, abs_isobel, abs_jsobel,
  178         [slice(None), slice(None), slice(1, None), slice(-1)],
  179         [slice(1, None), slice(-1), slice(1, None), slice(-1)],
  180         local_maxima)
  181     # ----- 90 to 135 degrees ------
  182     # Mix anti-diagonal and vertical
  183     #
  184     pts_plus = is_down & is_right
  185     pts_minus = is_up & is_left
  186     pts = ((pts_plus | pts_minus) & is_vertical)
  187     local_maxima = _set_local_maxima(
  188         magnitude, pts, abs_isobel, abs_jsobel,
  189         [slice(None), slice(None), slice(-1), slice(1, None)],
  190         [slice(1, None), slice(-1), slice(1, None), slice(-1)],
  191         local_maxima)
  192     # ----- 135 to 180 degrees ------
  193     # Mix anti-diagonal and anti-horizontal
  194     #
  195     pts = ((pts_plus | pts_minus) & is_horizontal)
  196     local_maxima = _set_local_maxima(
  197         magnitude, pts, abs_jsobel, abs_isobel,
  198         [slice(-1), slice(1, None), slice(-1), slice(1, None)],
  199         [slice(None), slice(None), slice(1, None), slice(-1)],
  200         local_maxima)
  201 
  202     return local_maxima
  203 
  204 
  205 def canny(image, sigma=1., low_threshold=None, high_threshold=None,
  206           mask=None, use_quantiles=False, *, mode='constant', cval=0.0):
  207     """Edge filter an image using the Canny algorithm.
  208 
  209     Parameters
  210     ----------
  211     image : 2D array
  212         Grayscale input image to detect edges on; can be of any dtype.
  213     sigma : float, optional
  214         Standard deviation of the Gaussian filter.
  215     low_threshold : float, optional
  216         Lower bound for hysteresis thresholding (linking edges).
  217         If None, low_threshold is set to 10% of dtype's max.
  218     high_threshold : float, optional
  219         Upper bound for hysteresis thresholding (linking edges).
  220         If None, high_threshold is set to 20% of dtype's max.
  221     mask : array, dtype=bool, optional
  222         Mask to limit the application of Canny to a certain area.
  223     use_quantiles : bool, optional
  224         If ``True`` then treat low_threshold and high_threshold as
  225         quantiles of the edge magnitude image, rather than absolute
  226         edge magnitude values. If ``True`` then the thresholds must be
  227         in the range [0, 1].
  228     mode : str, {'reflect', 'constant', 'nearest', 'mirror', 'wrap'}
  229         The ``mode`` parameter determines how the array borders are
  230         handled during Gaussian filtering, where ``cval`` is the value when
  231         mode is equal to 'constant'.
  232     cval : float, optional
  233         Value to fill past edges of input if `mode` is 'constant'.
  234 
  235     Returns
  236     -------
  237     output : 2D array (image)
  238         The binary edge map.
  239 
  240     See also
  241     --------
  242     skimage.sobel
  243 
  244     Notes
  245     -----
  246     The steps of the algorithm are as follows:
  247 
  248     * Smooth the image using a Gaussian with ``sigma`` width.
  249 
  250     * Apply the horizontal and vertical Sobel operators to get the gradients
  251       within the image. The edge strength is the norm of the gradient.
  252 
  253     * Thin potential edges to 1-pixel wide curves. First, find the normal
  254       to the edge at each point. This is done by looking at the
  255       signs and the relative magnitude of the X-Sobel and Y-Sobel
  256       to sort the points into 4 categories: horizontal, vertical,
  257       diagonal and antidiagonal. Then look in the normal and reverse
  258       directions to see if the values in either of those directions are
  259       greater than the point in question. Use interpolation to get a mix of
  260       points instead of picking the one that's the closest to the normal.
  261 
  262     * Perform a hysteresis thresholding: first label all points above the
  263       high threshold as edges. Then recursively label any point above the
  264       low threshold that is 8-connected to a labeled point as an edge.
  265 
  266     References
  267     ----------
  268     .. [1] Canny, J., A Computational Approach To Edge Detection, IEEE Trans.
  269            Pattern Analysis and Machine Intelligence, 8:679-714, 1986
  270            :DOI:`10.1109/TPAMI.1986.4767851`
  271     .. [2] William Green's Canny tutorial
  272            https://en.wikipedia.org/wiki/Canny_edge_detector
  273 
  274     Examples
  275     --------
  276     >>> from skimage import feature
  277     >>> rng = np.random.default_rng()
  278     >>> # Generate noisy image of a square
  279     >>> im = np.zeros((256, 256))
  280     >>> im[64:-64, 64:-64] = 1
  281     >>> im += 0.2 * rng.random(im.shape)
  282     >>> # First trial with the Canny filter, with the default smoothing
  283     >>> edges1 = feature.canny(im)
  284     >>> # Increase the smoothing for better results
  285     >>> edges2 = feature.canny(im, sigma=3)
  286 
  287     """
  288 
  289     # Regarding masks, any point touching a masked point will have a gradient
  290     # that is "infected" by the masked point, so it's enough to erode the
  291     # mask by one and then mask the output. We also mask out the border points
  292     # because who knows what lies beyond the edge of the image?
  293 
  294     check_nD(image, 2)
  295     dtype_max = dtype_limits(image, clip_negative=False)[1]
  296 
  297     if low_threshold is None:
  298         low_threshold = 0.1
  299     elif use_quantiles:
  300         if not(0.0 <= low_threshold <= 1.0):
  301             raise ValueError("Quantile thresholds must be between 0 and 1.")
  302     else:
  303         low_threshold /= dtype_max
  304 
  305     if high_threshold is None:
  306         high_threshold = 0.2
  307     elif use_quantiles:
  308         if not(0.0 <= high_threshold <= 1.0):
  309             raise ValueError("Quantile thresholds must be between 0 and 1.")
  310     else:
  311         high_threshold /= dtype_max
  312 
  313     if high_threshold < low_threshold:
  314         raise ValueError("low_threshold should be lower then high_threshold")
  315 
  316     # Image filtering
  317     smoothed, eroded_mask = _preprocess(image, mask, sigma, mode, cval)
  318 
  319     # Gradient magnitude estimation
  320     jsobel = ndi.sobel(smoothed, axis=1)
  321     isobel = ndi.sobel(smoothed, axis=0)
  322     magnitude = np.hypot(isobel, jsobel)
  323 
  324     if use_quantiles:
  325         low_threshold, high_threshold = np.percentile(magnitude,
  326                                                       [100.0 * low_threshold,
  327                                                        100.0 * high_threshold])
  328 
  329     # Non-maximum suppression
  330     local_maxima = _get_local_maxima(isobel, jsobel, magnitude, eroded_mask)
  331 
  332     # Double thresholding and edge traking
  333     low_mask = local_maxima & (magnitude >= low_threshold)
  334 
  335     #
  336     # Segment the low-mask, then only keep low-segments that have
  337     # some high_mask component in them
  338     #
  339     strel = np.ones((3, 3), bool)
  340     labels, count = ndi.label(low_mask, strel)
  341     if count == 0:
  342         return low_mask
  343 
  344     high_mask = local_maxima & (magnitude >= high_threshold)
  345     nonzero_sums = np.unique(labels[high_mask])
  346     good_label = np.zeros((count + 1,), bool)
  347     good_label[nonzero_sums] = True
  348     output_mask = good_label[labels]
  349     return output_mask