From 63b5c5a4a0bee2aea185b31af67da97fe751a2c2 Mon Sep 17 00:00:00 2001
From: "Josh Warner (Mac)" <warner.joshua@mayo.edu>
Date: Mon, 19 Nov 2012 23:38:58 -0600
Subject: [PATCH 01/11] 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 ------------------------------

From 226eaaff9250da4e73e4f0589cea41a13d961130 Mon Sep 17 00:00:00 2001
From: "Josh Warner (Mac)" <warner.joshua@mayo.edu>
Date: Tue, 20 Nov 2012 00:12:47 -0600
Subject: [PATCH 02/11] FEAT - Unified peak finder backend

Detailed changelist:
  * Fully documented new API for peak_local_max()
  * Removed commented-out old code for is_local_maximum()
  * Added "See also" and "Notes" to is_local_maximum() redirecting to peak_local_max()
---
 skimage/feature/peak.py         | 58 +++++++++++++++++++++-------
 skimage/morphology/watershed.py | 67 ++++++---------------------------
 2 files changed, 56 insertions(+), 69 deletions(-)

diff --git a/skimage/feature/peak.py b/skimage/feature/peak.py
index b1148c9a..cd740656 100644
--- a/skimage/feature/peak.py
+++ b/skimage/feature/peak.py
@@ -4,10 +4,10 @@ import scipy.ndimage as ndi
 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.
+    """
+    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`).
+    Peaks are the local maxima
 
     NOTE: If peaks are flat (i.e. multiple pixels have exact same intensity),
     the coordinates of all pixels are returned.
@@ -16,22 +16,54 @@ def peak_local_max(image, min_distance=10, threshold_abs=0, threshold_rel=0.1,
     ----------
     image : ndarray of floats
         Input image.
-    min_distance : int
-        Minimum number of pixels separating peaks and image boundary.
-    threshold : float
-        Deprecated. See `threshold_rel`.
-    threshold_abs : float
+
+    min_distance : int, default 10.
+        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 points, use `min_distance=1`.
+
+    threshold_abs : float, default 0.
         Minimum intensity of peaks.
-    threshold_rel : float
+
+    threshold_rel : float, default 0.1
         Minimum intensity of peaks calculated as `max(image) * threshold_rel`.
-    num_peaks : int
+
+    exclude_border : bool, default True
+        If True, `min_distance` excludes peaks from the border of the image as
+        well as from each other.
+
+    indices : bool, default True
+        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, default np.inf
         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` is 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.
+
+    threshold : float, optional
+        Deprecated. If provided as a kwarg, will override `threshold_rel`.
+        See `threshold_rel`.
 
     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
     -----
diff --git a/skimage/morphology/watershed.py b/skimage/morphology/watershed.py
index 840d4c90..f4c5a816 100644
--- a/skimage/morphology/watershed.py
+++ b/skimage/morphology/watershed.py
@@ -248,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))
@@ -281,63 +291,8 @@ 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,
                           footprint=footprint, indices=False,
                           exclude_border=False)

From 33c95140626374a68a8547114052aa25b8f3630d Mon Sep 17 00:00:00 2001
From: "Josh Warner (Mac)" <warner.joshua@mayo.edu>
Date: Tue, 20 Nov 2012 00:36:12 -0600
Subject: [PATCH 03/11] BUG - Fixed a bug where, when `indices` = False, output
 was not Boolean.

---
 skimage/feature/peak.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/skimage/feature/peak.py b/skimage/feature/peak.py
index cd740656..ffa942c4 100644
--- a/skimage/feature/peak.py
+++ b/skimage/feature/peak.py
@@ -119,7 +119,7 @@ def peak_local_max(image, min_distance=10, threshold_abs=0, threshold_rel=0.1,
         if indices is True:
             return np.transpose(out.nonzero())
         else:
-            return out
+            return out.astype(bool)
 
 
     if np.all(image == image.flat[0]):
@@ -164,6 +164,6 @@ def peak_local_max(image, min_distance=10, threshold_abs=0, threshold_rel=0.1,
     if indices is True:
         return coordinates
     else:
-        out = np.zeros_like(image)
-        out[coordinates[:, 0], coordinates[:, 1]] = 1
+        out = np.zeros_like(image, dtype=bool)
+        out[coordinates[:, 0], coordinates[:, 1]] = True
         return out

From 233b368195d6beaeba84794b5bfd8663e2bf45f2 Mon Sep 17 00:00:00 2001
From: "Josh Warner (Mac)" <warner.joshua@mayo.edu>
Date: Tue, 20 Nov 2012 00:45:01 -0600
Subject: [PATCH 04/11] BUG - Fixed a bug where the image was not properly
 masked for each unique label.

---
 skimage/feature/peak.py | 15 ++++++++-------
 1 file changed, 8 insertions(+), 7 deletions(-)

diff --git a/skimage/feature/peak.py b/skimage/feature/peak.py
index ffa942c4..fcd088c6 100644
--- a/skimage/feature/peak.py
+++ b/skimage/feature/peak.py
@@ -108,13 +108,14 @@ def peak_local_max(image, min_distance=10, threshold_abs=0, threshold_rel=0.1,
 
         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)
+            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,
+                                  **kwargs)
 
         if indices is True:
             return np.transpose(out.nonzero())

From da5a2a82213e01fc9b9d98f31cf740c6d3aa7767 Mon Sep 17 00:00:00 2001
From: "Josh Warner (Mac)" <warner.joshua@mayo.edu>
Date: Tue, 20 Nov 2012 01:11:31 -0600
Subject: [PATCH 05/11] BUG - Now parsing all unique labels != 0 properly, and
 set threshold_rel=0 in is_local_maximum() wrapper.

---
 skimage/feature/peak.py         | 4 +++-
 skimage/morphology/watershed.py | 4 ++--
 2 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/skimage/feature/peak.py b/skimage/feature/peak.py
index fcd088c6..ea2ace3e 100644
--- a/skimage/feature/peak.py
+++ b/skimage/feature/peak.py
@@ -106,8 +106,10 @@ def peak_local_max(image, min_distance=10, threshold_abs=0, threshold_rel=0.1,
             labels[mask] = rank_order(labels[mask])[0].astype(labels.dtype)
         labels = labels.astype(np.int32)
 
+        # New values for new ordering
+        label_values = np.unique(labels)
         out = np.zeros_like(image)
-        for label in 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,
diff --git a/skimage/morphology/watershed.py b/skimage/morphology/watershed.py
index f4c5a816..b6ae2d07 100644
--- a/skimage/morphology/watershed.py
+++ b/skimage/morphology/watershed.py
@@ -294,8 +294,8 @@ def is_local_maximum(image, labels=None, footprint=None):
 
     """
     return peak_local_max(image, labels=labels, min_distance=1,
-                          footprint=footprint, indices=False,
-                          exclude_border=False)
+                          threshold_rel=0, footprint=footprint,
+                          indices=False, exclude_border=False)
 
 
 # ---------------------- deprecated ------------------------------

From e7288f9140b1a1a8a545f19434ec5805006b4a1a Mon Sep 17 00:00:00 2001
From: "Josh Warner (Mac)" <warner.joshua@mayo.edu>
Date: Tue, 20 Nov 2012 01:17:42 -0600
Subject: [PATCH 06/11] BUG - Dtype conversion to Boolean in a new code path
 not included in current tests.

---
 skimage/feature/peak.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/skimage/feature/peak.py b/skimage/feature/peak.py
index ea2ace3e..2d40ded6 100644
--- a/skimage/feature/peak.py
+++ b/skimage/feature/peak.py
@@ -129,7 +129,7 @@ def peak_local_max(image, min_distance=10, threshold_abs=0, threshold_rel=0.1,
         if indices is True:
             return []
         else:
-            return np.zeros_like(image)
+            return np.zeros_like(image, dtype=bool)
 
     image = image.copy()
     # Non maximum filter

From d2262227d92bd135d44c712afa36f4660cf29425 Mon Sep 17 00:00:00 2001
From: "Josh Warner (Mac)" <warner.joshua@mayo.edu>
Date: Sat, 24 Nov 2012 17:24:14 -0600
Subject: [PATCH 07/11] BUG - `rank_order` now imported.

Also changed (peak.py):
  * Standardized documentation as requested.
  * Removed `threshold` as an optional kwarg.
  * Removed extra line break incorrect by PEP8 standards.

Also changed (watershed.py):
  * Added @deprecated decorator and import statement
---
 skimage/feature/peak.py         | 51 +++++++++++----------------------
 skimage/morphology/watershed.py |  2 ++
 2 files changed, 19 insertions(+), 34 deletions(-)

diff --git a/skimage/feature/peak.py b/skimage/feature/peak.py
index 2d40ded6..181dd600 100644
--- a/skimage/feature/peak.py
+++ b/skimage/feature/peak.py
@@ -1,9 +1,11 @@
 import numpy as np
 import scipy.ndimage as ndi
+from ..filter import rank_order
+
 
 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):
+                   footprint=None, labels=None):
     """
     Find peaks in an image, and return them as coordinates or a boolean array.
 
@@ -16,48 +18,34 @@ def peak_local_max(image, min_distance=10, threshold_abs=0, threshold_rel=0.1,
     ----------
     image : ndarray of floats
         Input image.
-
-    min_distance : int, default 10.
+    min_distance : int
         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.
+        `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 points, use `min_distance=1`.
-
-    threshold_abs : float, default 0.
+    threshold_abs : float
         Minimum intensity of peaks.
-
-    threshold_rel : float, default 0.1
+    threshold_rel : float
         Minimum intensity of peaks calculated as `max(image) * threshold_rel`.
-
-    exclude_border : bool, default True
+    exclude_border : bool
         If True, `min_distance` excludes peaks from the border of the image as
         well as from each other.
-
-    indices : bool, default True
+    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, default np.inf
+        with peaks present at True elements.
+    num_peaks : int
         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`.
-        Overrides `min_distance`, except for border exclusion if
-            `exclude_border` is True.
-
+        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.
 
-    threshold : float, optional
-        Deprecated. If provided as a kwarg, will override `threshold_rel`.
-        See `threshold_rel`.
-
     Returns
     -------
     output : (N, 2) array or ndarray of bools
@@ -116,15 +104,14 @@ def peak_local_max(image, min_distance=10, threshold_abs=0, threshold_rel=0.1,
                                   threshold_rel=threshold_rel,
                                   exclude_border=exclude_border,
                                   indices=False, num_peaks=np.inf,
-                                  footprint=footprint, labels=None,
-                                  **kwargs)
+                                  footprint=footprint, labels=None)
+            del maskim
 
         if indices is True:
             return np.transpose(out.nonzero())
         else:
             return out.astype(bool)
 
-
     if np.all(image == image.flat[0]):
         if indices is True:
             return []
@@ -149,15 +136,11 @@ def peak_local_max(image, min_distance=10, threshold_abs=0, threshold_rel=0.1,
         image[:, :min_distance] = 0
         image[:, -min_distance:] = 0
 
-    if kwargs.has_key('threshold'):
-        threshold_rel = kwargs['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]]
diff --git a/skimage/morphology/watershed.py b/skimage/morphology/watershed.py
index b6ae2d07..5c4ffd1a 100644
--- a/skimage/morphology/watershed.py
+++ b/skimage/morphology/watershed.py
@@ -29,6 +29,7 @@ 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
 
@@ -226,6 +227,7 @@ def watershed(image, markers, connectivity=None, offset=None, mask=None):
         return c_output
 
 
+@deprecated('filter.peak_local_max')
 def is_local_maximum(image, labels=None, footprint=None):
     """
     Return a boolean array of points that are local maxima

From 45a91aa9c5d9bdc2184440badd9ae2478fb296bd Mon Sep 17 00:00:00 2001
From: "Josh Warner (Mac)" <warner.joshua@mayo.edu>
Date: Sun, 25 Nov 2012 22:56:45 -0600
Subject: [PATCH 08/11] doc: deprecated string and docstring formatting fixes
 in `is_local_maximum`

---
 skimage/morphology/watershed.py | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/skimage/morphology/watershed.py b/skimage/morphology/watershed.py
index 5c4ffd1a..6d5a8371 100644
--- a/skimage/morphology/watershed.py
+++ b/skimage/morphology/watershed.py
@@ -227,7 +227,7 @@ def watershed(image, markers, connectivity=None, offset=None, mask=None):
         return c_output
 
 
-@deprecated('filter.peak_local_max')
+@deprecated('feature.peak_local_max')
 def is_local_maximum(image, labels=None, footprint=None):
     """
     Return a boolean array of points that are local maxima
@@ -236,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

From 98ce4751134b01f15c153a9f23ca6bdd48a259c3 Mon Sep 17 00:00:00 2001
From: "Josh Warner (Mac)" <warner.joshua@mayo.edu>
Date: Sun, 25 Nov 2012 23:09:06 -0600
Subject: [PATCH 09/11] style: initialize `out` at beginning of function

* remove API compatibility comment
* remove unnecessary `del maskim`
---
 skimage/feature/peak.py | 11 ++++-------
 1 file changed, 4 insertions(+), 7 deletions(-)

diff --git a/skimage/feature/peak.py b/skimage/feature/peak.py
index 181dd600..f68bd925 100644
--- a/skimage/feature/peak.py
+++ b/skimage/feature/peak.py
@@ -83,9 +83,9 @@ def peak_local_max(image, min_distance=10, threshold_abs=0, threshold_rel=0.1,
     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; for API compatibility with
-    # ..watershed.is_local_maximum()
+    # operating on each label separately
     if labels is not None:
         label_values = np.unique(labels)
         # Reorder label values to have consecutive integers (no gaps)
@@ -96,7 +96,6 @@ def peak_local_max(image, min_distance=10, threshold_abs=0, threshold_rel=0.1,
 
         # New values for new ordering
         label_values = np.unique(labels)
-        out = np.zeros_like(image)
         for label in label_values[label_values != 0]:
             maskim = (labels == label)
             out += peak_local_max(image * maskim, min_distance=min_distance,
@@ -105,18 +104,17 @@ def peak_local_max(image, min_distance=10, threshold_abs=0, threshold_rel=0.1,
                                   exclude_border=exclude_border,
                                   indices=False, num_peaks=np.inf,
                                   footprint=footprint, labels=None)
-            del maskim
 
         if indices is True:
             return np.transpose(out.nonzero())
         else:
-            return out.astype(bool)
+            return out.astype(np.bool)
 
     if np.all(image == image.flat[0]):
         if indices is True:
             return []
         else:
-            return np.zeros_like(image, dtype=bool)
+            return out
 
     image = image.copy()
     # Non maximum filter
@@ -150,6 +148,5 @@ def peak_local_max(image, min_distance=10, threshold_abs=0, threshold_rel=0.1,
     if indices is True:
         return coordinates
     else:
-        out = np.zeros_like(image, dtype=bool)
         out[coordinates[:, 0], coordinates[:, 1]] = True
         return out

From 4f68f1cd36ea3844fc7d8101e8416eca85379381 Mon Sep 17 00:00:00 2001
From: "Josh Warner (Mac)" <warner.joshua@mayo.edu>
Date: Mon, 3 Dec 2012 20:23:01 -0600
Subject: [PATCH 10/11] doc: fix and clarify `peak_local_max` docstring

---
 skimage/feature/peak.py | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/skimage/feature/peak.py b/skimage/feature/peak.py
index f68bd925..85d2cb4a 100644
--- a/skimage/feature/peak.py
+++ b/skimage/feature/peak.py
@@ -9,10 +9,11 @@ def peak_local_max(image, min_distance=10, threshold_abs=0, threshold_rel=0.1,
     """
     Find peaks in an image, and return them as coordinates or a boolean array.
 
-    Peaks are the local maxima
+    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
     ----------
@@ -23,7 +24,7 @@ def peak_local_max(image, min_distance=10, threshold_abs=0, threshold_rel=0.1,
         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 points, use `min_distance=1`.
+        To find the maximum number of peaks, use `min_distance=1`.
     threshold_abs : float
         Minimum intensity of peaks.
     threshold_rel : float

From 5cba69a4fad59b062910dc7dca1c53b489c2a0a8 Mon Sep 17 00:00:00 2001
From: "Josh Warner (Mac)" <warner.joshua@mayo.edu>
Date: Sat, 8 Dec 2012 21:12:26 -0600
Subject: [PATCH 11/11] test: add additional tests for 100% coverage of new
 `peak_local_max`

Also fix minor bug in label reordering.
---
 skimage/feature/peak.py            |  4 +--
 skimage/feature/tests/test_peak.py | 49 +++++++++++++++++++++++++++++-
 2 files changed, 50 insertions(+), 3 deletions(-)

diff --git a/skimage/feature/peak.py b/skimage/feature/peak.py
index 85d2cb4a..4ff70be1 100644
--- a/skimage/feature/peak.py
+++ b/skimage/feature/peak.py
@@ -91,8 +91,8 @@ def peak_local_max(image, min_distance=10, threshold_abs=0, threshold_rel=0.1,
         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)
+            mask = labels >= 1
+            labels[mask] = 1 + rank_order(labels[mask])[0].astype(labels.dtype)
         labels = labels.astype(np.int32)
 
         # New values for new ordering
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()