Merge pull request #21 from MrBago/peak_finding

RF - Make peak finding return zero peaks when there are no good peaks to...

dipy/reconst/peaks.py

@@ -93,23 +93,25 @@ def peak_directions(odf, sphere, relative_peak_threshold=.5,
                     min_separation_angle=25, minmax_norm=True):
     """Get the directions of odf peaks
 
+    Peaks are defined as points on the odf that are greater than at least one
+    neighbor and greater than or equal to all neighbors. Peaks are sorted in
+    descending order by their values then filtered based on their relative size
+    and spacing on the sphere. An odf may have 0 peaks, for example if the odf
+    is perfectly isotropic.
+
     Parameters
     ----------
     odf : 1d ndarray
         The odf function evaluated on the vertices of `sphere`
     sphere : Sphere
         The Sphere providing discrete directions for evaluation.
-    relative_peak_threshold : float
-        Only return peaks greater than ``relative_peak_threshold * m`` where m
-        is the largest peak.
-    min_separation_angle : float in [0, 90] The minimum distance between
-        directions. If two peaks are too close only the larger of the two is
-        returned.
-    minmax_norm : bool
-        If True min max normalization is applied internaly in the odf for the
-        peak finding. The returned ``values`` will be on the initial odf. This
-        parameter affects the ``relative_peak_threshold`` parameter because now
-        ``m`` will be 1.
+    relative_peak_threshold : float in [0., 1.]
+        Only peaks greater than ``min + relative_peak_threshold * scale`` are
+        kept, where ``min = max(0, odf.min())`` and
+        ``scale = odf.max() - min``.
+    min_separation_angle : float in [0, 90]
+        The minimum distance between directions. If two peaks are too close
+        only the larger of the two is returned.
 
     Returns
     -------
@@ -125,29 +127,28 @@ def peak_directions(odf, sphere, relative_peak_threshold=.5,
     If the odf has any negative values, they will be clipped to zeros.
 
     """
-
-    odf = np.clip(odf, 0, np.inf)
     values, indices = local_maxima(odf, sphere.edges)
 
     # If there is only one peak return
-    if len(indices) == 1:
+    n = len(values)
+    if n == 0 or (values[0] < 0.):
+        return np.zeros((0, 3)), np.zeros(0), np.zeros(0, dtype=int)
+    elif n == 1:
         return sphere.vertices[indices], values, indices
 
-    if minmax_norm:
-        odf_min = odf.min()
-        odf_max = values[0]
-        #because of the relative threshold this algorithm will
-        #give the same peaks as if we divide (values - odf_min) with
-        #(odf_max - odf_min) or not so here we skip the division
-        #to increase speed
-        values_norm = (values - odf_min)
-    else:
-        values_norm = values
+    odf_min = odf.min()
+    odf_min = odf_min if (odf_min >= 0.) else 0.
+    # because of the relative threshold this algorithm will give the same peaks
+    # as if we divide (values - odf_min) with (odf_max - odf_min) or not so
+    # here we skip the division to increase speed
+    values_norm = (values - odf_min)
 
+    # Remove small peaks
     n = search_descending(values_norm, relative_peak_threshold)
     indices = indices[:n]
     directions = sphere.vertices[indices]
 
+    # Remove peaks too close together
     directions, uniq = remove_similar_vertices(directions,
                                                min_separation_angle,
                                                return_index=True)

dipy/reconst/recspeed.pyx

@@ -204,17 +204,12 @@ def search_descending(cnp.ndarray[cnp.float_t, ndim=1, mode='c'] a,
         The greatest index such that ``all(a[:i] >= relative_threshold *
         a[0])``.
 
-    Note
-    ----
-    This function will never return 0, 1 is returned if ``a[0] <
-    relative_threshold * a[0]`` or if ``len(a) == 0``.
-
     """
     if a.shape[0] == 0:
-        return 1
+        return 0
 
     cdef:
-        size_t left = 1
+        size_t left = 0
         size_t right = a.shape[0]
         size_t mid
         double threshold = relative_threshould * a[0]
@@ -383,12 +378,6 @@ cdef long _compare_neighbors(double[:] odf, cnp.uint16_t[:, :] edges,
             wpeak_ptr[count] = i
             count += 1
 
-    # If count == 0, all values of odf are equal, and point 0 is returned as a
-    # peak to satisfy the requirement that peak_values[0] == max(odf).
-    if count == 0:
-        count = 1
-        wpeak_ptr[0] = 0
-
     return count
 
 

dipy/reconst/tests/test_peak_finding.py

@@ -134,7 +134,7 @@ def test_search_descending():
     npt.assert_equal(search_descending(a[:1], .5), 1)
 
     # Test very small array
-    npt.assert_equal(search_descending(a[:0], 1.), 1)
+    npt.assert_equal(search_descending(a[:0], 1.), 0)
 
 
 if __name__ == '__main__':

dipy/reconst/tests/test_peaks.py

@@ -355,16 +355,16 @@ def test_difference_with_minmax():
 
     assert_equal(len(values_2), 3)
 
-    _, values_3, _ = peak_directions(odf_gt, sphere, .30, 25.,
-                                     minmax_norm=False)
+    # Setting the smallest value of the odf to zero is like running
+    # peak_directions without the odf_min correction.
+    odf_gt[odf_gt.argmin()] = 0.
+    _, values_3, _ = peak_directions(odf_gt, sphere, .30, 25.,)
 
     assert_equal(len(values_3), 4)
 
     # we show here that to actually get that noisy peak out we need to
     # increase the peak threshold considerably
-    directions, values_4, indices = peak_directions(odf_gt, sphere,
-                                                    .60, 25.,
-                                                    minmax_norm=False)
+    directions, values_4, indices = peak_directions(odf_gt, sphere, .60, 25.,)
 
     assert_equal(len(values_4), 3)
     assert_almost_equal(values_1, values_4)
@@ -379,7 +379,9 @@ def test_degenerative_cases():
     directions, values, indices = peak_directions(odf, sphere, .5, 25)
     print(directions, values, indices)
 
-    assert_equal(values[0], 0)
+    assert_equal(len(values), 0)
+    assert_equal(len(directions), 0)
+    assert_equal(len(indices), 0)
 
     odf = np.zeros(sphere.vertices.shape[0])
     odf[0] = 0.020
@@ -394,7 +396,7 @@ def test_degenerative_cases():
     directions, values, indices = peak_directions(odf, sphere, .5, 25)
     print(directions, values, indices)
 
-    assert_equal(values[0], 0)
+    assert_equal(len(values), 0)
 
     odf = np.zeros(sphere.vertices.shape[0])
     odf[0] = 0.020
@@ -406,21 +408,11 @@ def test_degenerative_cases():
 
     odf = np.ones(sphere.vertices.shape[0])
     odf += 0.1 * np.random.rand(odf.shape[0])
-
-    directions, values, indices = peak_directions(odf, sphere, .5, 25)
-    assert_equal(len(directions) > 1, True)
-
-    odf = np.ones(sphere.vertices.shape[0])
-    odf -= np.finfo(np.float).eps * np.random.rand(odf.shape[0])
-
     directions, values, indices = peak_directions(odf, sphere, .5, 25)
-    assert_equal(len(directions) > 10, True)
+    assert_true(all(values > values[0] * .5))
+    assert_array_equal(values, odf[indices])
 
     odf = np.ones(sphere.vertices.shape[0])
-    directions, values, indices = peak_directions(odf, sphere, .5, 25)
-    assert_equal(values[0], 1)
-    assert_equal(len(values), 1)
-
     odf[1:] = np.finfo(np.float).eps * np.random.rand(odf.shape[0] - 1)
     directions, values, indices = peak_directions(odf, sphere, .5, 25)
 

dipy/tracking/markov.py

@@ -226,6 +226,8 @@ def _closest_peak(peak_directions, prev_step, cos_similarity):
     """
     if prev_step is None:
         return peak_directions
+    if len(peak_directions) == 0:
+        return None
 
     peak_dots = np.dot(peak_directions, prev_step)
     closest_peak = abs(peak_dots).argmax()