diff --git a/skimage/feature/peak.py b/skimage/feature/peak.py index 4765974e..4ff70be1 100644 --- a/skimage/feature/peak.py +++ b/skimage/feature/peak.py @@ -1,46 +1,66 @@ -import warnings import numpy as np -from scipy import ndimage +import scipy.ndimage as ndi +from ..filter import rank_order -def peak_local_max(image, min_distance=10, threshold='deprecated', - threshold_abs=0, threshold_rel=0.1, num_peaks=np.inf): - """Return coordinates of peaks in an image. +def peak_local_max(image, min_distance=10, threshold_abs=0, threshold_rel=0.1, + exclude_border=True, indices=True, num_peaks=np.inf, + footprint=None, labels=None): + """ + Find peaks in an image, and return them as coordinates or a boolean array. Peaks are the local maxima in a region of `2 * min_distance + 1` (i.e. peaks are separated by at least `min_distance`). - NOTE: If peaks are flat (i.e. multiple pixels have exact same intensity), - the coordinates of all pixels are returned. + NOTE: If peaks are flat (i.e. multiple adjacent pixels have identical + intensities), the coordinates of all such pixels are returned. Parameters ---------- image : ndarray of floats Input image. min_distance : int - Minimum number of pixels separating peaks and image boundary. - threshold : float - Deprecated. See `threshold_rel`. + Minimum number of pixels separating peaks in a region of `2 * + min_distance + 1` (i.e. peaks are separated by at least + `min_distance`). If `exclude_border` is True, this value also excludes + a border `min_distance` from the image boundary. + To find the maximum number of peaks, use `min_distance=1`. threshold_abs : float Minimum intensity of peaks. threshold_rel : float Minimum intensity of peaks calculated as `max(image) * threshold_rel`. + exclude_border : bool + If True, `min_distance` excludes peaks from the border of the image as + well as from each other. + indices : bool + If True, the output will be a matrix representing peak coordinates. + If False, the output will be a boolean matrix shaped as `image.shape` + with peaks present at True elements. num_peaks : int Maximum number of peaks. When the number of peaks exceeds `num_peaks`, - return `num_peaks` coordinates based on peak intensity. + 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`. Overrides + `min_distance`, except for border exclusion if `exclude_border=True`. + 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. Returns ------- - coordinates : (N, 2) array - (row, column) coordinates of peaks. + output : (N, 2) array or ndarray of bools + If `exclude_border = True` : (row, column) coordinates of peaks. + If `exclude_border = 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 a image. A maximum filter is used for finding local maxima. This operation - dilates the original image. After comparison between dilated and original image, - peak_local_max function returns the coordinates of peaks where - dilated image = original. + The peak local maximum function returns the coordinates of local peaks + (maxima) in a image. A maximum filter is used for finding local maxima. + This operation dilates the original image. After comparison between + dilated and original image, peak_local_max function returns the + coordinates of peaks where dilated image = original. Examples -------- @@ -64,35 +84,70 @@ def peak_local_max(image, min_distance=10, threshold='deprecated', array([[3, 2]]) """ + out = np.zeros_like(image, dtype=np.bool) + # In the case of labels, recursively build and return an output + # operating on each label separately + if labels is not None: + label_values = np.unique(labels) + # Reorder label values to have consecutive integers (no gaps) + if np.any(np.diff(label_values) != 1): + mask = labels >= 1 + labels[mask] = 1 + rank_order(labels[mask])[0].astype(labels.dtype) + labels = labels.astype(np.int32) + + # New values for new ordering + label_values = np.unique(labels) + for label in label_values[label_values != 0]: + maskim = (labels == label) + out += peak_local_max(image * maskim, min_distance=min_distance, + threshold_abs=threshold_abs, + threshold_rel=threshold_rel, + exclude_border=exclude_border, + indices=False, num_peaks=np.inf, + footprint=footprint, labels=None) + + if indices is True: + return np.transpose(out.nonzero()) + else: + return out.astype(np.bool) + if np.all(image == image.flat[0]): - return [] + if indices is True: + return [] + else: + return out + image = image.copy() # Non maximum filter - size = 2 * min_distance + 1 - image_max = ndimage.maximum_filter(image, size=size, mode='constant') + if footprint is not None: + image_max = ndi.maximum_filter(image, footprint=footprint, + mode='constant') + else: + size = 2 * min_distance + 1 + image_max = ndi.maximum_filter(image, size=size, mode='constant') mask = (image == image_max) image *= mask - # Remove the image borders - image[:min_distance] = 0 - image[-min_distance:] = 0 - image[:, :min_distance] = 0 - image[:, -min_distance:] = 0 + if exclude_border: + # Remove the image borders + image[:min_distance] = 0 + image[-min_distance:] = 0 + image[:, :min_distance] = 0 + image[:, -min_distance:] = 0 - if not threshold == 'deprecated': - msg = "`threshold` parameter deprecated; use `threshold_rel instead." - warnings.warn(msg, DeprecationWarning) - threshold_rel = threshold # find top peak candidates above a threshold peak_threshold = max(np.max(image.ravel()) * threshold_rel, threshold_abs) - image_t = (image > peak_threshold) * 1 # get coordinates of peaks - coordinates = np.transpose(image_t.nonzero()) + coordinates = np.transpose((image > peak_threshold).nonzero()) if coordinates.shape[0] > num_peaks: intensities = image[coordinates[:, 0], coordinates[:, 1]] idx_maxsort = np.argsort(intensities)[::-1] coordinates = coordinates[idx_maxsort][:num_peaks] - return coordinates + if indices is True: + return coordinates + else: + out[coordinates[:, 0], coordinates[:, 1]] = True + return out diff --git a/skimage/feature/tests/test_peak.py b/skimage/feature/tests/test_peak.py index 13457781..3ef1f12d 100644 --- a/skimage/feature/tests/test_peak.py +++ b/skimage/feature/tests/test_peak.py @@ -1,9 +1,17 @@ import numpy as np from numpy.testing import assert_array_almost_equal as assert_close - +import scipy.ndimage from skimage.feature import peak +def test_trivial_case(): + trivial = np.zeros((25, 25)) + peak_indices = peak.peak_local_max(trivial, min_distance=1, indices=True) + assert not peak_indices # inherent boolean-ness of empty list + peaks = peak.peak_local_max(trivial, min_distance=1, indices=False) + assert (peaks.astype(np.bool) == trivial).all() + + def test_noisy_peaks(): peak_locations = [(7, 7), (7, 13), (13, 7), (13, 13)] @@ -70,6 +78,45 @@ def test_num_peaks(): assert (3, 5) in peaks_limited +def test_reorder_labels(): + np.random.seed(21) + image = np.random.uniform(size=(40, 60)) + i, j = np.mgrid[0:40, 0:60] + labels = 1 + (i >= 20) + (j >= 30) * 2 + labels[labels == 4] = 5 + i, j = np.mgrid[-3:4, -3:4] + footprint = (i * i + j * j <= 9) + expected = np.zeros(image.shape, float) + for imin, imax in ((0, 20), (20, 40)): + for jmin, jmax in ((0, 30), (30, 60)): + expected[imin:imax, jmin:jmax] = scipy.ndimage.maximum_filter( + image[imin:imax, jmin:jmax], footprint=footprint) + expected = (expected == image) + result = peak.peak_local_max(image, labels=labels, min_distance=1, + threshold_rel=0, footprint=footprint, + indices=False, exclude_border=False) + assert (result == expected).all() + + +def test_indices_with_labels(): + np.random.seed(21) + image = np.random.uniform(size=(40, 60)) + i, j = np.mgrid[0:40, 0:60] + labels = 1 + (i >= 20) + (j >= 30) * 2 + i, j = np.mgrid[-3:4, -3:4] + footprint = (i * i + j * j <= 9) + expected = np.zeros(image.shape, float) + for imin, imax in ((0, 20), (20, 40)): + for jmin, jmax in ((0, 30), (30, 60)): + expected[imin:imax, jmin:jmax] = scipy.ndimage.maximum_filter( + image[imin:imax, jmin:jmax], footprint=footprint) + expected = (expected == image) + result = peak.peak_local_max(image, labels=labels, min_distance=1, + threshold_rel=0, footprint=footprint, + indices=True, exclude_border=False) + assert (result == np.transpose(expected.nonzero())).all() + + if __name__ == '__main__': from numpy import testing testing.run_module_suite() diff --git a/skimage/morphology/watershed.py b/skimage/morphology/watershed.py index 8a08fafc..6d5a8371 100644 --- a/skimage/morphology/watershed.py +++ b/skimage/morphology/watershed.py @@ -28,6 +28,8 @@ from _heapq import heappush, heappop import numpy as np import scipy.ndimage from ..filter import rank_order +from ..feature import peak_local_max +from .._shared.utils import deprecated from . import _watershed @@ -225,6 +227,7 @@ def watershed(image, markers, connectivity=None, offset=None, mask=None): return c_output +@deprecated('feature.peak_local_max') def is_local_maximum(image, labels=None, footprint=None): """ Return a boolean array of points that are local maxima @@ -233,10 +236,8 @@ def is_local_maximum(image, labels=None, footprint=None): ---------- image: ndarray (2-D, 3-D, ...) intensity image - labels: ndarray, optional find maxima only within labels. Zero is reserved for background. - footprint: ndarray of bools, optional binary mask indicating the neighborhood to be examined `footprint` must be a matrix with odd dimensions, the center is taken @@ -247,6 +248,16 @@ def is_local_maximum(image, labels=None, footprint=None): result: ndarray of bools mask that is True for pixels that are local maxima of `image` + See also + -------- + skimage.feature.peak_local_max: Unified peak finding backend. + The more capable backend for finding local maxima. + + Notes + ----- + This function is now a wrapper for skimage.feature.peak_local_max() and is + retained only for convenience and backward compatibility. + Examples -------- >>> image = np.zeros((4, 4)) @@ -280,63 +291,11 @@ def is_local_maximum(image, labels=None, footprint=None): [False, True, False, True], [False, False, False, False], [False, True, False, True]], dtype=bool) + """ - if labels is None: - labels = np.ones(image.shape, dtype=np.uint8) - if footprint is None: - footprint = np.ones([3] * image.ndim, dtype=np.uint8) - assert((np.all(footprint.shape) & 1) == 1) - footprint = (footprint != 0) - footprint_extent = (np.array(footprint.shape) - 1) // 2 - if np.all(footprint_extent == 0): - return labels > 0 - result = (labels > 0).copy() - # - # Create a labels matrix with zeros at the borders that might be - # hit by the footprint. - # - big_labels = np.zeros(np.array(labels.shape) + footprint_extent * 2, - labels.dtype) - big_labels[[slice(fe, -fe) for fe in footprint_extent]] = labels - # - # Find the relative indexes of each footprint element - # - image_strides = np.array(image.strides) // image.dtype.itemsize - big_strides = np.array(big_labels.strides) // big_labels.dtype.itemsize - result_strides = np.array(result.strides) // result.dtype.itemsize - footprint_offsets = np.mgrid[[slice(-fe, fe + 1) for fe in footprint_extent]] - - fp_image_offsets = np.sum(image_strides[:, np.newaxis] * - footprint_offsets[:, footprint], 0) - fp_big_offsets = np.sum(big_strides[:, np.newaxis] * - footprint_offsets[:, footprint], 0) - # - # Get the index of each labeled pixel in the image and big_labels arrays - # - indexes = np.mgrid[[slice(0, x) for x in labels.shape]][:, labels > 0] - image_indexes = np.sum(image_strides[:, np.newaxis] * indexes, 0) - big_indexes = np.sum(big_strides[:, np.newaxis] * - (indexes + footprint_extent[:, np.newaxis]), 0) - result_indexes = np.sum(result_strides[:, np.newaxis] * indexes, 0) - # - # Now operate on the raveled images - # - big_labels_raveled = big_labels.ravel() - image_raveled = image.ravel() - result_raveled = result.ravel() - # - # A hit is a hit if the label at the offset matches the label at the pixel - # and if the intensity at the pixel is greater or equal to the intensity - # at the offset. - # - for fp_image_offset, fp_big_offset in zip(fp_image_offsets, fp_big_offsets): - same_label = (big_labels_raveled[big_indexes + fp_big_offset] == - big_labels_raveled[big_indexes]) - less_than = (image_raveled[image_indexes[same_label]] < - image_raveled[image_indexes[same_label] + fp_image_offset]) - result_raveled[result_indexes[same_label][less_than]] = False - - return result + return peak_local_max(image, labels=labels, min_distance=1, + threshold_rel=0, footprint=footprint, + indices=False, exclude_border=False) # ---------------------- deprecated ------------------------------