RecoBundles and SLR workflows

bin/dipy_labelsbundles

@@ -0,0 +1,9 @@
+#!python
+
+from __future__ import division, print_function
+
+from dipy.workflows.flow_runner import run_flow
+from dipy.workflows.segment import LabelsBundlesFlow
+
+if __name__ == "__main__":
+    run_flow(LabelsBundlesFlow())

bin/dipy_recobundles

@@ -0,0 +1,9 @@
+#!python
+
+from __future__ import division, print_function
+
+from dipy.workflows.flow_runner import run_flow
+from dipy.workflows.segment import RecoBundlesFlow
+
+if __name__ == "__main__":
+    run_flow(RecoBundlesFlow())

bin/dipy_slr

@@ -0,0 +1,9 @@
+#!python
+
+from __future__ import division, print_function
+
+from dipy.workflows.flow_runner import run_flow
+from dipy.workflows.align import SlrWithQbxFlow
+
+if __name__ == "__main__":
+    run_flow(SlrWithQbxFlow())

dipy/align/streamlinear.py

@@ -10,13 +10,15 @@
                                       center_streamlines,
                                       set_number_of_points,
                                       select_random_set_of_streamlines,
-                                      length)
-from dipy.segment.clustering import QuickBundles
+                                      length,
+                                      Streamlines)
+from dipy.segment.clustering import QuickBundles, qbx_and_merge
 from dipy.core.geometry import (compose_transformations,
                                 compose_matrix,
                                 decompose_matrix)
 from dipy.utils.six import string_types
 from time import time
+from dipy.tracking.streamline import Streamlines
 
 MAX_DIST = 1e10
 LOG_MAX_DIST = np.log(MAX_DIST)
@@ -693,8 +695,15 @@ def bundle_min_distance_asymmetric_fast(t, static, moving, block_size):
 
 
 def remove_clusters_by_size(clusters, min_size=0):
+
     by_size = lambda c: len(c) >= min_size
-    return filter(by_size, clusters)
+    ob = filter(by_size, clusters)
+
+    centroids = Streamlines()
+    for cluster in ob:
+        centroids.append(cluster.centroid)
+
+    return centroids
 
 
 def progressive_slr(static, moving, metric, x0, bounds,
@@ -828,18 +837,17 @@ def progressive_slr(static, moving, metric, x0, bounds,
 
     return slm
 
-
-def slr_with_qb(static, moving,
-                x0='affine',
-                rm_small_clusters=50,
-                maxiter=100,
-                select_random=None,
-                verbose=False,
-                greater_than=50,
-                less_than=250,
-                qb_thr=15,
-                nb_pts=20,
-                progressive=True, num_threads=None):
+def slr_with_qbx(static, moving,
+                 x0='affine',
+                 rm_small_clusters=50,
+                 maxiter=100,
+                 select_random=None,
+                 verbose=False,
+                 greater_than=50,
+                 less_than=250,
+                 qbx_thr=[40, 30, 20, 15],
+                 nb_pts=20,
+                 progressive=True, rng=None, num_threads=None):
     """ Utility function for registering large tractograms.
 
     For efficiency we apply the registration on cluster centroids and remove
@@ -865,6 +873,9 @@ def slr_with_qb(static, moving,
     options : None or dict,
         Extra options to be used with the selected method.
 
+    rng : RandomState
+        If None creates RandomState in function.
+
     num_threads : int
         Number of threads. If None (default) then all available threads
         will be used. Only metrics using OpenMP will use this variable.
@@ -886,6 +897,9 @@ def slr_with_qb(static, moving,
             bundles using local and global streamline-based registration and
             clustering, Neuroimage, 2017.
     """
+    if rng is None:
+        rng = np.random.RandomState()
+
     if verbose:
         print('Static streamlines size {}'.format(len(static)))
         print('Moving streamlines size {}'.format(len(moving)))
@@ -898,8 +912,8 @@ def check_range(streamline, gt=greater_than, lt=less_than):
             return False
 
     # TODO change this to the new Streamlines API
-    streamlines1 = [s for s in static if check_range(s)]
-    streamlines2 = [s for s in moving if check_range(s)]
+    streamlines1 = Streamlines(static[np.array([check_range(s) for s in static])])
+    streamlines2 = Streamlines(moving[np.array([check_range(s) for s in moving])])
 
     if verbose:
 
@@ -910,29 +924,39 @@ def check_range(streamline, gt=greater_than, lt=less_than):
 
     if select_random is not None:
         rstreamlines1 = select_random_set_of_streamlines(streamlines1,
-                                                         select_random)
+                                                         select_random,
+                                                         rng=rng)
     else:
         rstreamlines1 = streamlines1
 
     rstreamlines1 = set_number_of_points(rstreamlines1, nb_pts)
-    qb1 = QuickBundles(threshold=qb_thr)
-    rstreamlines1 = [s.astype('f4') for s in rstreamlines1]
-    cluster_map1 = qb1.cluster(rstreamlines1)
+
+    # qb1 = QuickBundles(threshold=qb_thr)
+    rstreamlines1._data.astype('f4')
+    '''#rstreamlines1 = [s.astype('f4') for s in rstreamlines1]
+    # cluster_map1 = qb1.cluster(rstreamlines1)'''
+
+    cluster_map1 = qbx_and_merge(rstreamlines1, thresholds=qbx_thr, rng=rng)
     clusters1 = remove_clusters_by_size(cluster_map1, rm_small_clusters)
-    qb_centroids1 = [cluster.centroid for cluster in clusters1]
+
+    qb_centroids1 = clusters1
 
     if select_random is not None:
         rstreamlines2 = select_random_set_of_streamlines(streamlines2,
-                                                         select_random)
+                                                         select_random,
+                                                         rng=rng)
     else:
         rstreamlines2 = streamlines2
 
     rstreamlines2 = set_number_of_points(rstreamlines2, nb_pts)
-    qb2 = QuickBundles(threshold=qb_thr)
-    rstreamlines2 = [s.astype('f4') for s in rstreamlines2]
-    cluster_map2 = qb2.cluster(rstreamlines2)
+    rstreamlines2._data.astype('f4')
+    '''# qb2 = QuickBundles(threshold=qb_thr)
+    #rstreamlines2 = [s.astype('f4') for s in rstreamlines2]
+    # cluster_map2 = qb2.cluster(rstreamlines2)'''
+    cluster_map2 = qbx_and_merge(rstreamlines2, thresholds=qbx_thr, rng=rng)
+
     clusters2 = remove_clusters_by_size(cluster_map2, rm_small_clusters)
-    qb_centroids2 = [cluster.centroid for cluster in clusters2]
+    qb_centroids2 = clusters2
 
     if verbose:
         t = time()
@@ -967,7 +991,7 @@ def check_range(streamline, gt=greater_than, lt=less_than):
 # Garyfallidis et al. Recognition of white matter
 # bundles using local and global streamline-based registration and
 # clustering, Neuroimage, 2017.
-whole_brain_slr = slr_with_qb
+whole_brain_slr = slr_with_qbx
 
 
 def _threshold(x, th):

dipy/align/tests/test_whole_brain_slr.py

@@ -4,7 +4,7 @@
                            assert_array_almost_equal)
 from dipy.data import get_data
 from dipy.tracking.streamline import Streamlines
-from dipy.align.streamlinear import whole_brain_slr, slr_with_qb
+from dipy.align.streamlinear import whole_brain_slr, slr_with_qbx
 from dipy.tracking.distances import bundles_distances_mam
 from dipy.align.streamlinear import transform_streamlines
 from dipy.align.streamlinear import compose_matrix44, decompose_matrix44
@@ -22,8 +22,9 @@ def test_whole_brain_slr():
     f2._data += np.array([50, 0, 0])
 
     moved, transform, qb_centroids1, qb_centroids2 = whole_brain_slr(
-            f1, f2, verbose=True, rm_small_clusters=2, greater_than=0,
-            less_than=np.inf, qb_thr=5, progressive=False)
+            f1, f2, x0='affine', verbose=True, rm_small_clusters=2,
+            greater_than=0, less_than=np.inf,
+            qbx_thr=[5, 2, 1], progressive=False)
 
     # we can check the quality of registration by comparing the matrices
     # MAM streamline distances before and after SLR
@@ -33,31 +34,37 @@ def test_whole_brain_slr():
     d12_minsum = np.sum(np.min(D12, axis=0))
     d1m_minsum = np.sum(np.min(D1M, axis=0))
 
+    print("distances= ", d12_minsum, " ", d1m_minsum)
+
     assert_equal(d1m_minsum < d12_minsum, True)
 
-    assert_array_almost_equal(transform[:3, 3], [-50, -0, -0], 3)
+    assert_array_almost_equal(transform[:3, 3], [-50, -0, -0], 2)
 
     # check rotation
+
     mat = compose_matrix44([0, 0, 0, 15, 0, 0])
 
     f3 = f.copy()
     f3 = transform_streamlines(f3, mat)
 
-    moved, transform, qb_centroids1, qb_centroids2 = slr_with_qb(
+    moved, transform, qb_centroids1, qb_centroids2 = slr_with_qbx(
             f1, f3, verbose=False, rm_small_clusters=1, greater_than=20,
-            less_than=np.inf, qb_thr=2, progressive=True)
+            less_than=np.inf, qbx_thr=[2],
+            progressive=True)
 
     # we can also check the quality by looking at the decomposed transform
+
     assert_array_almost_equal(decompose_matrix44(transform)[3], -15, 2)
 
-    moved, transform, qb_centroids1, qb_centroids2 = slr_with_qb(
+    moved, transform, qb_centroids1, qb_centroids2 = slr_with_qbx(
             f1, f3, verbose=False, rm_small_clusters=1, select_random=400,
-            greater_than=20,
-            less_than=np.inf, qb_thr=2, progressive=True)
+            greater_than=20, less_than=np.inf, qbx_thr=[2],
+            progressive=True)
 
     # we can also check the quality by looking at the decomposed transform
-    assert_array_almost_equal(decompose_matrix44(transform)[3], -15, 2)
 
+    assert_array_almost_equal(decompose_matrix44(transform)[3], -15, 2)
 
 if __name__ == '__main__':
-    run_module_suite()
+    # run_module_suite()
+    test_whole_brain_slr()

dipy/data/__init__.py

@@ -46,7 +46,11 @@
                                read_tissue_data,
                                fetch_cfin_multib,
                                read_cfin_dwi,
-                               read_cfin_t1)
+                               read_cfin_t1,
+                               fetch_target_tractogram_hcp,
+                               fetch_bundle_atlas_hcp842,
+                               read_bundle_atlas_hcp842,
+                               read_target_tractogram_hcp)
 
 from ..utils.arrfuncs import as_native_array
 from dipy.tracking.streamline import relist_streamlines

dipy/data/fetcher.py

@@ -1,5 +1,4 @@
 from __future__ import division, print_function, absolute_import
-
 import os
 import sys
 import contextlib
@@ -32,6 +31,7 @@
 UW_RW_URL = \
   "https://digital.lib.washington.edu/researchworks/bitstream/handle/"
 
+
 class FetcherError(Exception):
     pass
 
@@ -430,6 +430,28 @@ def fetcher():
          " More details about the data are available in their paper: " +
          " https://www.nature.com/articles/sdata201672"))
 
+fetch_bundle_atlas_hcp842 = _make_fetcher(
+    "fetch_bundle_atlas_hcp842",
+    pjoin(dipy_home, 'bundle_atlas_hcp842'),
+    'https://ndownloader.figshare.com/files/',
+    ['11921522'],
+    ['Atlas_in_MNI_Space_16_bundles.zip'],
+    ['b071f3e851f21ba1749c02fc6beb3118'],
+    doc="Download atlas tractogram from the hcp842 dataset with its bundles",
+    data_size="200MB",
+    unzip=True)
+
+fetch_target_tractogram_hcp = _make_fetcher(
+    "fetch_target_tractogram_hcp",
+    pjoin(dipy_home, 'target_tractogram_hcp'),
+    'https://ndownloader.figshare.com/files/',
+    ['12871127'],
+    ['hcp_tractogram.zip'],
+    ['fa25ef19c9d3748929b6423397963b6a'],
+    doc="Download tractogram of one of the hcp dataset subjects",
+    data_size="541MB",
+    unzip=True)
+
 
 def read_scil_b0():
     """ Load GE 3T b0 image form the scil b0 dataset.
@@ -1031,4 +1053,61 @@ def read_cfin_t1():
     """
     files, folder = fetch_cfin_multib()
     img = nib.load(pjoin(folder, 'T1.nii'))
-    return img, gtab
+    return img  # , gtab
+
+
+def read_bundle_atlas_hcp842():
+    """
+    Returns
+    -------
+    file1 : string
+    file2 : string
+    """
+    file1 = pjoin(dipy_home,
+                  'bundle_atlas_hcp842',
+                  'Atlas_in_MNI_Space_16_bundles',
+                  'whole_brain',
+                  'whole_brain_MNI.trk')
+
+    file2 = pjoin(dipy_home,
+                  'bundle_atlas_hcp842',
+                  'Atlas_in_MNI_Space_16_bundles',
+                  'bundles',
+                  '*.trk')
+
+    return file1, file2
+
+
+def read_two_hcp842_bundle():
+    """
+    Returns
+    -------
+    file1 : string
+    file2 : string
+    """
+    file1 = pjoin(dipy_home,
+                  'bundle_atlas_hcp842',
+                  'Atlas_in_MNI_Space_16_bundles',
+                  'bundles',
+                  'AF_L.trk')
+
+    file2 = pjoin(dipy_home,
+                  'bundle_atlas_hcp842',
+                  'Atlas_in_MNI_Space_16_bundles',
+                  'bundles',
+                  'CST_L.trk')
+
+    return file1, file2
+
+
+def read_target_tractogram_hcp():
+    """
+    Returns
+    -------
+    file1 : string
+    """
+    file1 = pjoin(dipy_home, 'target_tractogram_hcp',
+                              'hcp_tractogram',
+                              'streamlines.trk')
+
+    return file1

dipy/reconst/csdeconv.py

@@ -811,9 +811,9 @@ def auto_response(gtab, data, roi_center=None, roi_radius=10, fa_thr=0.7,
     fa_thr : float
         FA threshold
     fa_callable : callable
-        A callable that defines an operation that compares FA with the fa_thr. The operator
-        should have two positional arguments (e.g., `fa_operator(FA, fa_thr)`) and it should
-        return a bool array.
+        A callable that defines an operation that compares FA with the fa_thr.
+        The operator should have two positional arguments
+        (e.g., `fa_operator(FA, fa_thr)`) and it should return a bool array.
     return_number_of_voxels : bool
         If True, returns the number of voxels used for estimating the response
         function.