diff --git a/dipy/reconst/peaks.py b/dipy/reconst/peaks.py index 0699e8bc9f..e5151189e4 100644 --- a/dipy/reconst/peaks.py +++ b/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) diff --git a/dipy/reconst/recspeed.pyx b/dipy/reconst/recspeed.pyx index e70df2d2c4..3c3c4f14b8 100644 --- a/dipy/reconst/recspeed.pyx +++ b/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 diff --git a/dipy/reconst/tests/test_peak_finding.py b/dipy/reconst/tests/test_peak_finding.py index 7b827bf5a4..20c5c5ce7e 100644 --- a/dipy/reconst/tests/test_peak_finding.py +++ b/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__': diff --git a/dipy/reconst/tests/test_peaks.py b/dipy/reconst/tests/test_peaks.py index ce94eb24cb..c5f00d16f0 100644 --- a/dipy/reconst/tests/test_peaks.py +++ b/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) diff --git a/dipy/tracking/markov.py b/dipy/tracking/markov.py index 43ebc6adf5..db860f0a0a 100644 --- a/dipy/tracking/markov.py +++ b/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()