From 63b5c5a4a0bee2aea185b31af67da97fe751a2c2 Mon Sep 17 00:00:00 2001 From: "Josh Warner (Mac)" Date: Mon, 19 Nov 2012 23:38:58 -0600 Subject: [PATCH] FEAT - combined API from is_local_maximum() into peak_local_max() is_local_maximum() is a wrapper function for peak_local_max() is_local_maximum() runs much faster (~20% of prior runtime, nearly = to peak_local_max()) All tests in .feature and .morphology subpackages pass as written with these changes. Todo: * Fully document API * remove commented-out old algorithm in is_local_maximum() * add new tests for full coverage of new, more complex peak_local_max() --- skimage/feature/peak.py | 85 +++++++++++++++++------- skimage/morphology/watershed.py | 112 +++++++++++++++++--------------- 2 files changed, 120 insertions(+), 77 deletions(-) diff --git a/skimage/feature/peak.py b/skimage/feature/peak.py index 4765974e..b1148c9a 100644 --- a/skimage/feature/peak.py +++ b/skimage/feature/peak.py @@ -1,10 +1,9 @@ -import warnings import numpy as np -from scipy import ndimage +import scipy.ndimage as ndi - -def peak_local_max(image, min_distance=10, threshold='deprecated', - threshold_abs=0, threshold_rel=0.1, num_peaks=np.inf): +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, **kwargs): """Return coordinates of peaks in an image. Peaks are the local maxima in a region of `2 * min_distance + 1` @@ -36,11 +35,11 @@ def peak_local_max(image, min_distance=10, threshold='deprecated', 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,25 +63,60 @@ def peak_local_max(image, min_distance=10, threshold='deprecated', array([[3, 2]]) """ + # In the case of labels, recursively build and return an output + # operating on each label separately; for API compatibility with + # ..watershed.is_local_maximum() + 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 >= 0 + labels[mask] = rank_order(labels[mask])[0].astype(labels.dtype) + labels = labels.astype(np.int32) + + out = np.zeros_like(image) + for label in labels: + out += peak_local_max(image, 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, + **kwargs) + + if indices is True: + return np.transpose(out.nonzero()) + else: + return out + + if np.all(image == image.flat[0]): - return [] + if indices is True: + return [] + else: + return np.zeros_like(image) + 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 kwargs.has_key('threshold'): + threshold_rel = kwargs['threshold'] - 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 @@ -95,4 +129,9 @@ def peak_local_max(image, min_distance=10, threshold='deprecated', idx_maxsort = np.argsort(intensities)[::-1] coordinates = coordinates[idx_maxsort][:num_peaks] - return coordinates + if indices is True: + return coordinates + else: + out = np.zeros_like(image) + out[coordinates[:, 0], coordinates[:, 1]] = 1 + return out diff --git a/skimage/morphology/watershed.py b/skimage/morphology/watershed.py index 8a08fafc..840d4c90 100644 --- a/skimage/morphology/watershed.py +++ b/skimage/morphology/watershed.py @@ -28,6 +28,7 @@ from _heapq import heappush, heappop import numpy as np import scipy.ndimage from ..filter import rank_order +from ..feature import peak_local_max from . import _watershed @@ -281,62 +282,65 @@ def is_local_maximum(image, labels=None, footprint=None): [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]] + # 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 + # 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 result + return peak_local_max(image, labels=labels, min_distance=1, + footprint=footprint, indices=False, + exclude_border=False) # ---------------------- deprecated ------------------------------