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()