| 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420 |
- from warnings import warn
- import numpy as np
- import scipy.ndimage as ndi
- from .. import measure
- 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, kind="stable")
- coord = np.transpose(coord)[idx_maxsort]
- if np.isfinite(num_peaks):
- max_out = int(num_peaks)
- else:
- max_out = None
- if min_distance > 1:
- coord = ensure_spacing(
- coord, spacing=min_distance, p_norm=p_norm, max_out=max_out
- )
- 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='nearest')
- 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
- def peak_local_max(
- image,
- min_distance=1,
- threshold_abs=None,
- threshold_rel=None,
- exclude_border=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.
- 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 or None, optional
- Minimum intensity of peaks. By default, the absolute threshold is
- the minimum intensity of the image.
- threshold_rel : float or None, 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.
- 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
- The coordinates of the peaks.
- 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 images, this function returns the coordinates
- 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].copy()
- # 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
- )
- return coordinates
- 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.intensity_max)[::-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 neighborhood 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 neighborhood suppression
- ycoords_nh = ycoords_idx + ycoords_ext
- xcoords_nh = xcoords_idx + xcoords_ext
- # no reflection for distance neighborhood
- 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 neighborhood
- 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
|