mirror of
https://github.com/wassname/scikit-image.git
synced 2026-07-29 11:26:57 +08:00
Merge pull request #376 from JDWarner/unify_peak_finder_backend
ENH: Unify peak finder backend.
This commit is contained in:
Stefan van der Walt
+88
-33
@@ -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
|
||||
|
||||
@@ -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()
|
||||
|
||||
@@ -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 ------------------------------
|
||||
|
||||
Reference in New Issue
Block a user