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

@@ -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_abs=0, threshold_rel=0.1,
-def peak_local_max(image, min_distance=10, threshold='deprecated',
+                   exclude_border=True, indices=True, num_peaks=np.inf,
-                   threshold_abs=0, threshold_rel=0.1, 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)
+    The peak local maximum function returns the coordinates of local peaks
-    in a image. A maximum filter is used for finding local maxima. This operation
+    (maxima) in a image. A maximum filter is used for finding local maxima.
-    dilates the original image. After comparison between dilated and original image,
+    This operation dilates the original image. After comparison between
-    peak_local_max function returns the coordinates of peaks where
+    dilated and original image, peak_local_max function returns the
-    dilated image = original.
+    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
+    if footprint is not None:
-    image_max = ndimage.maximum_filter(image, size=size, mode='constant')
+        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
+    if exclude_border:
-    image[:min_distance] = 0
+        # Remove the image borders
-    image[-min_distance:] = 0
+        image[:min_distance] = 0
-    image[:, :min_distance] = 0
+        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

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:
+    # if labels is None:
-        labels = np.ones(image.shape, dtype=np.uint8)
+    #     labels = np.ones(image.shape, dtype=np.uint8)
-    if footprint is None:
+    # if footprint is None:
-        footprint = np.ones([3] * image.ndim, dtype=np.uint8)
+    #     footprint = np.ones([3] * image.ndim, dtype=np.uint8)
-    assert((np.all(footprint.shape) & 1) == 1)
+    # assert((np.all(footprint.shape) & 1) == 1)
-    footprint = (footprint != 0)
+    # footprint = (footprint != 0)
-    footprint_extent = (np.array(footprint.shape) - 1) // 2
+    # footprint_extent = (np.array(footprint.shape) - 1) // 2
-    if np.all(footprint_extent == 0):
+    # if np.all(footprint_extent == 0):
-        return labels > 0
+    #     return labels > 0
-    result = (labels > 0).copy()
+    # result = (labels > 0).copy()
-    #
+    # #
-    # Create a labels matrix with zeros at the borders that might be
+    # # Create a labels matrix with zeros at the borders that might be
-    # hit by the footprint.
+    # # hit by the footprint.
-    #
+    # #
-    big_labels = np.zeros(np.array(labels.shape) + footprint_extent * 2,
+    # big_labels = np.zeros(np.array(labels.shape) + footprint_extent * 2,
-                          labels.dtype)
+    #                       labels.dtype)
-    big_labels[[slice(fe, -fe) for fe in footprint_extent]] = labels
+    # big_labels[[slice(fe, -fe) for fe in footprint_extent]] = labels
-    #
+    # #
-    # Find the relative indexes of each footprint element
+    # # Find the relative indexes of each footprint element
-    #
+    # #
-    image_strides = np.array(image.strides) // image.dtype.itemsize
+    # image_strides = np.array(image.strides) // image.dtype.itemsize
-    big_strides = np.array(big_labels.strides) // big_labels.dtype.itemsize
+    # big_strides = np.array(big_labels.strides) // big_labels.dtype.itemsize
-    result_strides = np.array(result.strides) // result.dtype.itemsize
+    # result_strides = np.array(result.strides) // result.dtype.itemsize
-    footprint_offsets = np.mgrid[[slice(-fe, fe + 1) for fe in footprint_extent]]
+    # footprint_offsets = np.mgrid[[slice(-fe, fe + 1) for fe in footprint_extent]]
 
-    fp_image_offsets = np.sum(image_strides[:, np.newaxis] *
+    # fp_image_offsets = np.sum(image_strides[:, np.newaxis] *
-                              footprint_offsets[:, footprint], 0)
+    #                           footprint_offsets[:, footprint], 0)
-    fp_big_offsets = np.sum(big_strides[:, np.newaxis] *
+    # fp_big_offsets = np.sum(big_strides[:, np.newaxis] *
-                            footprint_offsets[:, footprint], 0)
+    #                         footprint_offsets[:, footprint], 0)
-    #
+    # #
-    # Get the index of each labeled pixel in the image and big_labels arrays
+    # # 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]
+    # indexes = np.mgrid[[slice(0, x) for x in labels.shape]][:, labels > 0]
-    image_indexes = np.sum(image_strides[:, np.newaxis] * indexes, 0)
+    # image_indexes = np.sum(image_strides[:, np.newaxis] * indexes, 0)
-    big_indexes = np.sum(big_strides[:, np.newaxis] *
+    # big_indexes = np.sum(big_strides[:, np.newaxis] *
-                         (indexes + footprint_extent[:, np.newaxis]), 0)
+    #                      (indexes + footprint_extent[:, np.newaxis]), 0)
-    result_indexes = np.sum(result_strides[:, np.newaxis] * indexes, 0)
+    # result_indexes = np.sum(result_strides[:, np.newaxis] * indexes, 0)
-    #
+    # #
-    # Now operate on the raveled images
+    # # Now operate on the raveled images
-    #
+    # #
-    big_labels_raveled = big_labels.ravel()
+    # big_labels_raveled = big_labels.ravel()
-    image_raveled = image.ravel()
+    # image_raveled = image.ravel()
-    result_raveled = result.ravel()
+    # result_raveled = result.ravel()
-    #
+    # #
-    # A hit is a hit if the label at the offset matches the label at the pixel
+    # # 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
+    # # and if the intensity at the pixel is greater or equal to the intensity
-    # at the offset.
+    # # at the offset.
-    #
+    # #
-    for fp_image_offset, fp_big_offset in zip(fp_image_offsets, fp_big_offsets):
+    # 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] ==
+    #     same_label = (big_labels_raveled[big_indexes + fp_big_offset] ==
-                      big_labels_raveled[big_indexes])
+    #                   big_labels_raveled[big_indexes])
-        less_than = (image_raveled[image_indexes[same_label]] <
+    #     less_than = (image_raveled[image_indexes[same_label]] <
-                     image_raveled[image_indexes[same_label] + fp_image_offset])
+    #                  image_raveled[image_indexes[same_label] + fp_image_offset])
-        result_raveled[result_indexes[same_label][less_than]] = False
+    #     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 ------------------------------