diff --git a/skimage/feature/peak.py b/skimage/feature/peak.py
index 225f42cd..280263af 100644
--- a/skimage/feature/peak.py
+++ b/skimage/feature/peak.py
@@ -1,29 +1,46 @@
+import warnings
 import numpy as np
 from scipy import ndimage
 
 
-def peak_local_max(image, min_distance=10, threshold=0.1):
+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.
 
     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.
+
     Parameters
     ----------
     image: ndarray of floats
         Input image.
 
-    min_distance: int, optional
+    min_distance: int
         Minimum number of pixels separating peaks and image boundary.
 
-    threshold: float, optional
-        Candidate peaks are calculated as `max(image) * threshold`.
+    threshold : float
+        Deprecated. See `threshold_rel`.
+
+    threshold_abs: float
+        Minimum intensity of peaks.
+
+    threshold_rel: float
+        Minimum intensity of peaks calculated as `max(image) * threshold_rel`.
+
+    num_peaks : int
+        Maximum number of peaks. When the number of peaks exceeds `num_peaks`,
+        return `num_peaks` coordinates based on peak intensity.
 
     Returns
     -------
     coordinates : (N, 2) array
         (row, column) coordinates of peaks.
     """
+    if np.all(image == image.flat[0]):
+        return []
     image = image.copy()
     # Non maximum filter
     size = 2 * min_distance + 1
@@ -37,12 +54,21 @@ def peak_local_max(image, min_distance=10, threshold=0.1):
     image[:, :min_distance] = 0
     image[:, -min_distance:] = 0
 
-    # find top corner candidates above a threshold
-    corner_threshold = np.max(image.ravel()) * threshold
-    image_t = (image >= corner_threshold) * 1
+    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())
 
+    if len(coordinates) > num_peaks:
+        intensities = image[tuple(coordinates.T)]
+        idx_maxsort = np.argsort(intensities)[::-1]
+        coordinates = coordinates[idx_maxsort][:2]
+
     return coordinates
 
diff --git a/skimage/feature/tests/test_peak.py b/skimage/feature/tests/test_peak.py
index 19f9e95f..91d38d99 100644
--- a/skimage/feature/tests/test_peak.py
+++ b/skimage/feature/tests/test_peak.py
@@ -1,6 +1,7 @@
 import numpy as np
+from numpy.testing import assert_array_almost_equal as assert_close
 
-from skimage import feature
+from skimage.feature import peak
 
 
 def test_noisy_peaks():
@@ -11,13 +12,56 @@ def test_noisy_peaks():
     for r, c in peak_locations:
         image[r, c] = 1
 
-    peaks_detected = feature.peak_local_max(image, min_distance=5)
+    peaks_detected = peak.peak_local_max(image, min_distance=5)
 
     assert len(peaks_detected) == len(peak_locations)
     for loc in peaks_detected:
         assert tuple(loc) in peak_locations
 
 
+def test_relative_threshold():
+    image = np.zeros((5, 5), dtype=np.uint8)
+    image[1, 1] = 10
+    image[3, 3] = 20
+    peaks = peak.peak_local_max(image, min_distance=1, threshold_rel=0.5)
+    assert len(peaks) == 1
+    assert_close(peaks, [(3, 3)])
+
+
+def test_absolute_threshold():
+    image = np.zeros((5, 5), dtype=np.uint8)
+    image[1, 1] = 10
+    image[3, 3] = 20
+    peaks = peak.peak_local_max(image, min_distance=1, threshold_abs=10)
+    assert len(peaks) == 1
+    assert_close(peaks, [(3, 3)])
+
+
+def test_constant_image():
+    image = 128 * np.ones((20, 20), dtype=np.uint8)
+    peaks = peak.peak_local_max(image, min_distance=1)
+    assert len(peaks) == 0
+
+
+def test_flat_peak():
+    image = np.zeros((5, 5), dtype=np.uint8)
+    image[1:3, 1:3] = 10
+    peaks = peak.peak_local_max(image, min_distance=1)
+    assert len(peaks) == 4
+
+
+def test_num_peaks():
+    image = np.zeros((3, 7), dtype=np.uint8)
+    image[1, 1] = 10
+    image[1, 3] = 11
+    image[1, 5] = 12
+    assert len(peak.peak_local_max(image, min_distance=1)) == 3
+    peaks_limited = peak.peak_local_max(image, min_distance=1, num_peaks=2)
+    assert len(peaks_limited) == 2
+    assert (1, 3) in peaks_limited
+    assert (1, 5) in peaks_limited
+
+
 if __name__ == '__main__':
     from numpy import testing
     testing.run_module_suite()