from warnings import warn import numpy as np import scipy.ndimage as ndi from .. import measure from .._shared.utils import remove_arg from .._shared.coord import ensure_spacing def _get_high_intensity_peaks(image, mask, num_peaks, min_distance, p_norm): """ Return the highest intensity peak coordinates. """ # get coordinates of peaks coord = np.nonzero(mask) intensities = image[coord] # Highest peak first idx_maxsort = np.argsort(-intensities) coord = np.transpose(coord)[idx_maxsort] coord = ensure_spacing(coord, spacing=min_distance, p_norm=p_norm) if len(coord) > num_peaks: coord = coord[:num_peaks] return coord def _get_peak_mask(image, footprint, threshold, mask=None): """ Return the mask containing all peak candidates above thresholds. """ if footprint.size == 1 or image.size == 1: return image > threshold image_max = ndi.maximum_filter(image, footprint=footprint, mode='constant') out = image == image_max # no peak for a trivial image image_is_trivial = np.all(out) if mask is None else np.all(out[mask]) if image_is_trivial: out[:] = False if mask is not None: # isolated pixels in masked area are returned as peaks isolated_px = np.logical_xor(mask, ndi.binary_opening(mask)) out[isolated_px] = True out &= image > threshold return out def _exclude_border(label, border_width): """Set label border values to 0. """ # zero out label borders for i, width in enumerate(border_width): if width == 0: continue label[(slice(None),) * i + (slice(None, width),)] = 0 label[(slice(None),) * i + (slice(-width, None),)] = 0 return label def _get_threshold(image, threshold_abs, threshold_rel): """Return the threshold value according to an absolute and a relative value. """ threshold = threshold_abs if threshold_abs is not None else image.min() if threshold_rel is not None: threshold = max(threshold, threshold_rel * image.max()) return threshold def _get_excluded_border_width(image, min_distance, exclude_border): """Return border_width values relative to a min_distance if requested. """ if isinstance(exclude_border, bool): border_width = (min_distance if exclude_border else 0,) * image.ndim elif isinstance(exclude_border, int): if exclude_border < 0: raise ValueError("`exclude_border` cannot be a negative value") border_width = (exclude_border,) * image.ndim elif isinstance(exclude_border, tuple): if len(exclude_border) != image.ndim: raise ValueError( "`exclude_border` should have the same length as the " "dimensionality of the image.") for exclude in exclude_border: if not isinstance(exclude, int): raise ValueError( "`exclude_border`, when expressed as a tuple, must only " "contain ints." ) if exclude < 0: raise ValueError( "`exclude_border` can not be a negative value") border_width = exclude_border else: raise TypeError( "`exclude_border` must be bool, int, or tuple with the same " "length as the dimensionality of the image.") return border_width @remove_arg("indices", changed_version="0.20") def peak_local_max(image, min_distance=1, threshold_abs=None, threshold_rel=None, exclude_border=True, indices=True, num_peaks=np.inf, footprint=None, labels=None, num_peaks_per_label=np.inf, p_norm=np.inf): """Find peaks in an image as coordinate list or boolean mask. Peaks are the local maxima in a region of `2 * min_distance + 1` (i.e. peaks are separated by at least `min_distance`). If both `threshold_abs` and `threshold_rel` are provided, the maximum of the two is chosen as the minimum intensity threshold of peaks. .. versionchanged:: 0.18 Prior to version 0.18, peaks of the same height within a radius of `min_distance` were all returned, but this could cause unexpected behaviour. From 0.18 onwards, an arbitrary peak within the region is returned. See issue gh-2592. Parameters ---------- image : ndarray Input image. min_distance : int, optional The minimal allowed distance separating peaks. To find the maximum number of peaks, use `min_distance=1`. threshold_abs : float, optional Minimum intensity of peaks. By default, the absolute threshold is the minimum intensity of the image. threshold_rel : float, optional Minimum intensity of peaks, calculated as `max(image) * threshold_rel`. exclude_border : int, tuple of ints, or bool, optional If positive integer, `exclude_border` excludes peaks from within `exclude_border`-pixels of the border of the image. If tuple of non-negative ints, the length of the tuple must match the input array's dimensionality. Each element of the tuple will exclude peaks from within `exclude_border`-pixels of the border of the image along that dimension. If True, takes the `min_distance` parameter as value. If zero or False, peaks are identified regardless of their distance from the border. indices : bool, optional If True, the output will be an array representing peak coordinates. The coordinates are sorted according to peaks values (Larger first). If False, the output will be a boolean array shaped as `image.shape` with peaks present at True elements. ``indices`` is deprecated and will be removed in version 0.20. Default behavior will be to always return peak coordinates. You can obtain a mask as shown in the example below. num_peaks : int, optional Maximum number of peaks. When the number of peaks exceeds `num_peaks`, return `num_peaks` peaks based on highest peak intensity. footprint : ndarray of bools, optional If provided, `footprint == 1` represents the local region within which to search for peaks at every point in `image`. labels : ndarray of ints, optional If provided, each unique region `labels == value` represents a unique region to search for peaks. Zero is reserved for background. num_peaks_per_label : int, optional Maximum number of peaks for each label. p_norm : float Which Minkowski p-norm to use. Should be in the range [1, inf]. A finite large p may cause a ValueError if overflow can occur. ``inf`` corresponds to the Chebyshev distance and 2 to the Euclidean distance. Returns ------- output : ndarray or ndarray of bools * If `indices = True` : (row, column, ...) coordinates of peaks. * If `indices = False` : Boolean array shaped like `image`, with peaks represented by True values. Notes ----- The peak local maximum function returns the coordinates of local peaks (maxima) in an image. Internally, a maximum filter is used for finding local maxima. This operation dilates the original image. After comparison of the dilated and original image, this function returns the coordinates or a mask of the peaks where the dilated image equals the original image. See also -------- skimage.feature.corner_peaks Examples -------- >>> img1 = np.zeros((7, 7)) >>> img1[3, 4] = 1 >>> img1[3, 2] = 1.5 >>> img1 array([[0. , 0. , 0. , 0. , 0. , 0. , 0. ], [0. , 0. , 0. , 0. , 0. , 0. , 0. ], [0. , 0. , 0. , 0. , 0. , 0. , 0. ], [0. , 0. , 1.5, 0. , 1. , 0. , 0. ], [0. , 0. , 0. , 0. , 0. , 0. , 0. ], [0. , 0. , 0. , 0. , 0. , 0. , 0. ], [0. , 0. , 0. , 0. , 0. , 0. , 0. ]]) >>> peak_local_max(img1, min_distance=1) array([[3, 2], [3, 4]]) >>> peak_local_max(img1, min_distance=2) array([[3, 2]]) >>> img2 = np.zeros((20, 20, 20)) >>> img2[10, 10, 10] = 1 >>> img2[15, 15, 15] = 1 >>> peak_idx = peak_local_max(img2, exclude_border=0) >>> peak_idx array([[10, 10, 10], [15, 15, 15]]) >>> peak_mask = np.zeros_like(img2, dtype=bool) >>> peak_mask[tuple(peak_idx.T)] = True >>> np.argwhere(peak_mask) array([[10, 10, 10], [15, 15, 15]]) """ if (footprint is None or footprint.size == 1) and min_distance < 1: warn("When min_distance < 1, peak_local_max acts as finding " "image > max(threshold_abs, threshold_rel * max(image)).", RuntimeWarning, stacklevel=2) border_width = _get_excluded_border_width(image, min_distance, exclude_border) threshold = _get_threshold(image, threshold_abs, threshold_rel) if footprint is None: size = 2 * min_distance + 1 footprint = np.ones((size, ) * image.ndim, dtype=bool) else: footprint = np.asarray(footprint) if labels is None: # Non maximum filter mask = _get_peak_mask(image, footprint, threshold) mask = _exclude_border(mask, border_width) # Select highest intensities (num_peaks) coordinates = _get_high_intensity_peaks(image, mask, num_peaks, min_distance, p_norm) else: _labels = _exclude_border(labels.astype(int, casting="safe"), border_width) if np.issubdtype(image.dtype, np.floating): bg_val = np.finfo(image.dtype).min else: bg_val = np.iinfo(image.dtype).min # For each label, extract a smaller image enclosing the object of # interest, identify num_peaks_per_label peaks labels_peak_coord = [] for label_idx, roi in enumerate(ndi.find_objects(_labels)): if roi is None: continue # Get roi mask label_mask = labels[roi] == label_idx + 1 # Extract image roi img_object = image[roi] # Ensure masked values don't affect roi's local peaks img_object[np.logical_not(label_mask)] = bg_val mask = _get_peak_mask(img_object, footprint, threshold, label_mask) coordinates = _get_high_intensity_peaks(img_object, mask, num_peaks_per_label, min_distance, p_norm) # transform coordinates in global image indices space for idx, s in enumerate(roi): coordinates[:, idx] += s.start labels_peak_coord.append(coordinates) if labels_peak_coord: coordinates = np.vstack(labels_peak_coord) else: coordinates = np.empty((0, 2), dtype=int) if len(coordinates) > num_peaks: out = np.zeros_like(image, dtype=bool) out[tuple(coordinates.T)] = True coordinates = _get_high_intensity_peaks(image, out, num_peaks, min_distance, p_norm) if indices: return coordinates else: out = np.zeros_like(image, dtype=bool) out[tuple(coordinates.T)] = True return out def _prominent_peaks(image, min_xdistance=1, min_ydistance=1, threshold=None, num_peaks=np.inf): """Return peaks with non-maximum suppression. Identifies most prominent features separated by certain distances. Non-maximum suppression with different sizes is applied separately in the first and second dimension of the image to identify peaks. Parameters ---------- image : (M, N) ndarray Input image. min_xdistance : int Minimum distance separating features in the x dimension. min_ydistance : int Minimum distance separating features in the y dimension. threshold : float Minimum intensity of peaks. Default is `0.5 * max(image)`. num_peaks : int Maximum number of peaks. When the number of peaks exceeds `num_peaks`, return `num_peaks` coordinates based on peak intensity. Returns ------- intensity, xcoords, ycoords : tuple of array Peak intensity values, x and y indices. """ img = image.copy() rows, cols = img.shape if threshold is None: threshold = 0.5 * np.max(img) ycoords_size = 2 * min_ydistance + 1 xcoords_size = 2 * min_xdistance + 1 img_max = ndi.maximum_filter1d(img, size=ycoords_size, axis=0, mode='constant', cval=0) img_max = ndi.maximum_filter1d(img_max, size=xcoords_size, axis=1, mode='constant', cval=0) mask = (img == img_max) img *= mask img_t = img > threshold label_img = measure.label(img_t) props = measure.regionprops(label_img, img_max) # Sort the list of peaks by intensity, not left-right, so larger peaks # in Hough space cannot be arbitrarily suppressed by smaller neighbors props = sorted(props, key=lambda x: x.max_intensity)[::-1] coords = np.array([np.round(p.centroid) for p in props], dtype=int) img_peaks = [] ycoords_peaks = [] xcoords_peaks = [] # relative coordinate grid for local neighbourhood suppression ycoords_ext, xcoords_ext = np.mgrid[-min_ydistance:min_ydistance + 1, -min_xdistance:min_xdistance + 1] for ycoords_idx, xcoords_idx in coords: accum = img_max[ycoords_idx, xcoords_idx] if accum > threshold: # absolute coordinate grid for local neighbourhood suppression ycoords_nh = ycoords_idx + ycoords_ext xcoords_nh = xcoords_idx + xcoords_ext # no reflection for distance neighbourhood ycoords_in = np.logical_and(ycoords_nh > 0, ycoords_nh < rows) ycoords_nh = ycoords_nh[ycoords_in] xcoords_nh = xcoords_nh[ycoords_in] # reflect xcoords and assume xcoords are continuous, # e.g. for angles: # (..., 88, 89, -90, -89, ..., 89, -90, -89, ...) xcoords_low = xcoords_nh < 0 ycoords_nh[xcoords_low] = rows - ycoords_nh[xcoords_low] xcoords_nh[xcoords_low] += cols xcoords_high = xcoords_nh >= cols ycoords_nh[xcoords_high] = rows - ycoords_nh[xcoords_high] xcoords_nh[xcoords_high] -= cols # suppress neighbourhood img_max[ycoords_nh, xcoords_nh] = 0 # add current feature to peaks img_peaks.append(accum) ycoords_peaks.append(ycoords_idx) xcoords_peaks.append(xcoords_idx) img_peaks = np.array(img_peaks) ycoords_peaks = np.array(ycoords_peaks) xcoords_peaks = np.array(xcoords_peaks) if num_peaks < len(img_peaks): idx_maxsort = np.argsort(img_peaks)[::-1][:num_peaks] img_peaks = img_peaks[idx_maxsort] ycoords_peaks = ycoords_peaks[idx_maxsort] xcoords_peaks = xcoords_peaks[idx_maxsort] return img_peaks, xcoords_peaks, ycoords_peaks