Merge pull request #2098 from sciunto/bug_localmaxpeak

FIX: select num_peaks if labels is specified

skimage/feature/peak.py

@@ -3,9 +3,26 @@
 from ..filters import rank_order
 
 
+def _get_high_intensity_peaks(image, mask, num_peaks):
+    """
+    Return the highest intensity peak coordinates.
+    """
+    # get coordinates of peaks
+    coord = np.nonzero(mask)
+    # select num_peaks peaks
+    if len(coord[0]) > num_peaks:
+        intensities = image[coord]
+        idx_maxsort = np.argsort(intensities)
+        coord = np.transpose(coord)[idx_maxsort][-num_peaks:]
+    else:
+        coord = np.column_stack(coord)
+    return coord
+
+
 def peak_local_max(image, min_distance=1, threshold_abs=None,
                    threshold_rel=None, exclude_border=True, indices=True,
-                   num_peaks=np.inf, footprint=None, labels=None):
+                   num_peaks=np.inf, footprint=None, labels=None,
+                   num_peaks_per_label=np.inf):
     """Find peaks in an image as coordinate list or boolean mask.
 
     Peaks are the local maxima in a region of `2 * min_distance + 1`
@@ -48,6 +65,8 @@ def peak_local_max(image, min_distance=1, threshold_abs=None,
     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.
+    num_peaks_per_label : int, optional
+        Maximum number of peaks for each label.
 
     Returns
     -------
@@ -92,7 +111,6 @@ def peak_local_max(image, min_distance=1, threshold_abs=None,
     array([[10, 10, 10]])
 
     """
-
     if type(exclude_border) == bool:
         exclude_border = min_distance if exclude_border else 0
 
@@ -116,13 +134,18 @@ def peak_local_max(image, min_distance=1, threshold_abs=None,
                                   threshold_abs=threshold_abs,
                                   threshold_rel=threshold_rel,
                                   exclude_border=exclude_border,
-                                  indices=False, num_peaks=np.inf,
+                                  indices=False, num_peaks=num_peaks_per_label,
                                   footprint=footprint, labels=None)
 
+        # Select highest intensities (num_peaks)
+        coordinates = _get_high_intensity_peaks(image, out, num_peaks)
+
         if indices is True:
-            return np.transpose(out.nonzero())
+            return coordinates
         else:
-            return out.astype(np.bool)
+            nd_indices = tuple(coordinates.T)
+            out[nd_indices] = True
+            return out
 
     if np.all(image == image.flat[0]):
         if indices is True:
@@ -158,14 +181,8 @@ def peak_local_max(image, min_distance=1, threshold_abs=None,
     if thresholds:
         mask &= image > max(thresholds)
 
-    # get coordinates of peaks
-    coordinates = np.transpose(mask.nonzero())
-
-    if coordinates.shape[0] > num_peaks:
-        intensities = image.flat[np.ravel_multi_index(coordinates.transpose(),
-                                                      image.shape)]
-        idx_maxsort = np.argsort(intensities)[::-1]
-        coordinates = coordinates[idx_maxsort][:num_peaks]
+    # Select highest intensities (num_peaks)
+    coordinates = _get_high_intensity_peaks(image, mask, num_peaks)
 
     if indices is True:
         return coordinates

skimage/feature/tests/test_peak.py

@@ -5,9 +5,6 @@
 from skimage.feature import peak
 
 
-np.random.seed(21)
-
-
 def test_trivial_case():
     trivial = np.zeros((25, 25))
     peak_indices = peak.peak_local_max(trivial, min_distance=1, indices=True)
@@ -21,6 +18,7 @@ def test_noisy_peaks():
     peak_locations = [(7, 7), (7, 13), (13, 7), (13, 13)]
 
     # image with noise of amplitude 0.8 and peaks of amplitude 1
+    np.random.seed(21)
     image = 0.8 * np.random.rand(20, 20)
     for r, c in peak_locations:
         image[r, c] = 1
@@ -85,6 +83,47 @@ def test_num_peaks():
     assert (3, 5) in peaks_limited
 
 
+def test_num_peaks_and_labels():
+    image = np.zeros((7, 7), dtype=np.uint8)
+    labels = np.zeros((7, 7), dtype=np.uint8) + 20
+    image[1, 1] = 10
+    image[1, 3] = 11
+    image[1, 5] = 12
+    image[3, 5] = 8
+    image[5, 3] = 7
+    peaks_limited = peak.peak_local_max(
+        image, min_distance=1, threshold_abs=0, labels=labels)
+    assert len(peaks_limited) == 5
+    peaks_limited = peak.peak_local_max(
+        image, min_distance=1, threshold_abs=0, labels=labels, num_peaks=2)
+    assert len(peaks_limited) == 2
+
+
+def test_num_peaks_tot_vs_labels_4quadrants():
+    np.random.seed(21)
+    image = np.random.uniform(size=(20, 30))
+    i, j = np.mgrid[0:20, 0:30]
+    labels = 1 + (i >= 10) + (j >= 15) * 2
+    result = peak.peak_local_max(image, labels=labels,
+                                 min_distance=1, threshold_rel=0,
+                                 indices=True,
+                                 num_peaks=np.inf,
+                                 num_peaks_per_label=2)
+    assert len(result) == 8
+    result = peak.peak_local_max(image, labels=labels,
+                                 min_distance=1, threshold_rel=0,
+                                 indices=True,
+                                 num_peaks=np.inf,
+                                 num_peaks_per_label=1)
+    assert len(result) == 4
+    result = peak.peak_local_max(image, labels=labels,
+                                 min_distance=1, threshold_rel=0,
+                                 indices=True,
+                                 num_peaks=2,
+                                 num_peaks_per_label=2)
+    assert len(result) == 2
+
+
 def test_num_peaks3D():
     # Issue 1354: the old code only hold for 2D arrays
     # and this code would die with IndexError
@@ -95,6 +134,7 @@ def test_num_peaks3D():
 
 
 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
@@ -114,6 +154,7 @@ def test_reorder_labels():
 
 
 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
@@ -272,14 +313,15 @@ def test_adjacent_different_objects():
 
 
 def test_four_quadrants():
-    image = np.random.uniform(size=(40, 60))
-    i, j = np.mgrid[0:40, 0:60]
-    labels = 1 + (i >= 20) + (j >= 30) * 2
+    np.random.seed(21)
+    image = np.random.uniform(size=(20, 30))
+    i, j = np.mgrid[0:20, 0:30]
+    labels = 1 + (i >= 10) + (j >= 15) * 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)):
+    for imin, imax in ((0, 10), (10, 20)):
+        for jmin, jmax in ((0, 15), (15, 30)):
             expected[imin:imax, jmin:jmax] = ndi.maximum_filter(
                 image[imin:imax, jmin:jmax], footprint=footprint)
     expected = (expected == image)
@@ -293,6 +335,7 @@ def test_disk():
     '''regression test of img-1194, footprint = [1]
     Test peak.peak_local_max when every point is a local maximum
     '''
+    np.random.seed(21)
     image = np.random.uniform(size=(10, 20))
     footprint = np.array([[1]])
     result = peak.peak_local_max(image, labels=np.ones((10, 20)),