Merge pull request #538 from bloomdt-uw/enh-callosum-rois

Enh Add Callosum ROIs support

AFQ/api.py

@@ -53,6 +53,10 @@ def do_preprocessing():
 BUNDLES = ["ATR", "CGC", "CST", "IFO", "ILF", "SLF", "ARC", "UNC",
            "FA", "FP"]
 
+CALLOSUM_BUNDLES = ["AntFrontal", "Motor", "Occipital", "Orbital",
+                    "PostParietal", "SupFrontal", "SupParietal",
+                    "Temporal"]
+
 # See: https://www.cmu.edu/dietrich/psychology/cognitiveaxon/documents/yeh_etal_2018.pdf  # noqa
 
 RECO_BUNDLES_16 = [
@@ -73,7 +77,9 @@ def do_preprocessing():
 DIPY_GH = "https://github.com/dipy/dipy/blob/master/dipy/"
 
 
-def make_bundle_dict(bundle_names=BUNDLES, seg_algo="afq", resample_to=False):
+def make_bundle_dict(bundle_names=BUNDLES,
+                     seg_algo="afq",
+                     resample_to=False):
     """
     Create a bundle dictionary, needed for the segmentation
 
@@ -94,7 +100,26 @@ def make_bundle_dict(bundle_names=BUNDLES, seg_algo="afq", resample_to=False):
         If set, templates will be resampled to the affine and shape of this
         image.
     """
+    logger = logging.getLogger('AFQ.api')
     if seg_algo == "afq":
+        if "FP" in bundle_names and "Occipital" in bundle_names:
+            logger.warning((
+                f"FP and Occipital bundles are co-located, and AFQ"
+                f" assigns each streamline to only one bundle."
+                f" Only Occipital will be used."))
+            bundle_names.remove("FP")
+        if "FA" in bundle_names and "Orbital" in bundle_names:
+            logger.warning((
+                f"FA and Orbital bundles are co-located, and AFQ"
+                f" assigns each streamline to only one bundle."
+                f" Only Orbital will be used."))
+            bundle_names.remove("FA")
+        if "FA" in bundle_names and "AntFrontal" in bundle_names:
+            logger.warning((
+                f"FA and AntFrontal bundles are co-located, and AFQ"
+                f" assigns each streamline to only one bundle."
+                f" Only AntFrontal will be used."))
+            bundle_names.remove("FA")
         templates = afd.read_templates(resample_to=resample_to)
         callosal_templates = afd.read_callosum_templates(
             resample_to=resample_to)
@@ -121,6 +146,15 @@ def make_bundle_dict(bundle_names=BUNDLES, seg_algo="afq", resample_to=False):
                     'cross_midline': True,
                     'uid': uid}
                 uid += 1
+            elif name in CALLOSUM_BUNDLES:
+                afq_bundles[name] = {
+                    'ROIs': [callosal_templates["L_" + name],
+                             callosal_templates["R_" + name],
+                             callosal_templates["Callosum_midsag"]],
+                    'rules': [True, True, True],
+                    'cross_midline': True,
+                    'uid': uid}
+                uid += 1
             # SLF is a special case, because it has an exclusion ROI:
             elif name == "SLF":
                 for hemi in ['_R', '_L']:
@@ -135,16 +169,20 @@ def make_bundle_dict(bundle_names=BUNDLES, seg_algo="afq", resample_to=False):
                     uid += 1
             else:
                 for hemi in ['_R', '_L']:
-                    afq_bundles[name + hemi] = {
-                        'ROIs': [templates[name + '_roi1' + hemi],
-                                 templates[name + '_roi2' + hemi]],
-                        'rules': [True, True],
-                        'prob_map': templates[name + hemi + '_prob_map'],
-                        'cross_midline': False,
-                        'uid': uid}
+                    if (templates.get(name + '_roi1' + hemi)
+                            and templates.get(name + '_roi2' + hemi)
+                            and templates.get(name + hemi + '_prob_map')):
+                        afq_bundles[name + hemi] = {
+                            'ROIs': [templates[name + '_roi1' + hemi],
+                                     templates[name + '_roi2' + hemi]],
+                            'rules': [True, True],
+                            'prob_map': templates[name + hemi + '_prob_map'],
+                            'cross_midline': False,
+                            'uid': uid}
+                    else:
+                        logger.warning(f"{name} is not in AFQ templates")
 
                     uid += 1
-
     elif seg_algo.startswith("reco"):
         if seg_algo.endswith("80"):
             bundle_dict = afd.read_hcp_atlas(80)

AFQ/data.py

@@ -1733,7 +1733,7 @@ def organize_cfin_data(path=None):
            "PipelineDescription": {"Name": "dipy"}})
 
 
-def organize_stanford_data(path=None):
+def organize_stanford_data(path=None, clear_previous_afq=False):
     """
     If necessary, downloads the Stanford HARDI dataset into DIPY directory and
     creates a BIDS compliant file-system structure in AFQ data directory:
@@ -1759,6 +1759,8 @@ def organize_stanford_data(path=None):
                         ├── sub-01_ses-01_dwi.bvec
                         └── sub-01_ses-01_dwi.nii.gz
 
+    If clear_previous_afq is True and there is an afq folder in derivatives,
+    it will be removed.
     """
     logger = logging.getLogger('AFQ.data')
 
@@ -1777,6 +1779,11 @@ def organize_stanford_data(path=None):
     dmriprep_folder = op.join(derivatives_path, 'vistasoft')
     freesurfer_folder = op.join(derivatives_path, 'freesurfer')
 
+    if clear_previous_afq:
+        afq_folder = op.join(derivatives_path, 'afq')
+        if op.exists(afq_folder):
+            shutil.rmtree(afq_folder)
+
     if not op.exists(derivatives_path):
         logger.info(f'creating derivatives directory: {derivatives_path}')
 
@@ -2092,7 +2099,15 @@ def bundles_to_aal(bundles, atlas=None):
         "UNC_L": [['leftanttemporal'], ['leftuncinatefront']],
         "UNC_R": [['rightanttemporal'], ['rightuncinatefront']],
         "ARC_L": [['leftfrontal'], ['leftarctemp']],
-        "ARC_R": [['rightfrontal'], ['rightarctemp']]}
+        "ARC_R": [['rightfrontal'], ['rightarctemp']],
+        "AntFrontal": [None, None],
+        "Motor": [None, None],
+        "Occipital": [None, None],
+        "Orbital": [None, None],
+        "PostParietal": [None, None],
+        "SupFrontal": [None, None],
+        "SupParietal": [None, None],
+        "Temporal": [None, None]}
 
     targets = []
     for bundle in bundles:

AFQ/tests/test_api.py

@@ -174,6 +174,47 @@ def create_dummy_bids_path(n_subjects, n_sessions, share_sessions=True):
     return bids_dir
 
 
+def test_make_bundle_dict():
+    """
+    Tests bundle dict
+    """
+
+    # test defaults
+    afq_bundles = api.make_bundle_dict()
+
+    # bundles restricted within hemisphere
+    # NOTE: FA and FP cross midline so are removed
+    # NOTE: all others generate two bundles
+    num_hemi_bundles = (len(api.BUNDLES)-2)*2
+
+    # bundles that cross the midline
+    num_whole_bundles = 2
+
+    assert len(afq_bundles) == num_hemi_bundles + num_whole_bundles
+
+    # Arcuate Fasciculus
+    afq_bundles = api.make_bundle_dict(bundle_names=["ARC"])
+
+    assert len(afq_bundles) == 2
+
+    # Forceps Minor
+    afq_bundles = api.make_bundle_dict(bundle_names=["FA"])
+
+    assert len(afq_bundles) == 1
+
+    # Cingulum Hippocampus
+    # not included but exists in templates
+    afq_bundles = api.make_bundle_dict(bundle_names=["HCC"])
+
+    assert len(afq_bundles) == 2
+
+    # Vertical Occipital Fasciculus
+    # not included and does not exist in templates
+    afq_bundles = api.make_bundle_dict(bundle_names=["VOF"])
+
+    assert len(afq_bundles) == 0
+
+
 @pytest.mark.nightly4
 def test_AFQ_custom_tract():
     """

examples/plot_afq_api.py

@@ -38,8 +38,12 @@
 #
 # This data represents the required preprocessed diffusion data necessary for
 # intializing the AFQ object (which we will do next)
+#
+# The clear_previous_afq is used to remove any previous runs of the afq object
+# stored in the AFQ_data/stanford_hardi/ BIDS directory. Set it to false if
+# you want to use the results of previous runs. 
 
-afd.organize_stanford_data()
+afd.organize_stanford_data(clear_previous_afq=True)
 
 ##########################################################################
 # Initialize an AFQ object:

examples/plot_afq_callosal.py

@@ -0,0 +1,60 @@
+"""
+==========================
+Callosal bundles using AFQ API
+==========================
+An example using the AFQ API to find callosal bundles using the templates from:
+http://hdl.handle.net/1773/34926
+"""
+import os.path as op
+
+import plotly
+
+from AFQ import api
+from AFQ.mask import RoiMask
+import AFQ.data as afd
+
+##########################################################################
+# Get some example data
+# ---------------------
+#
+# Retrieves `Stanford HARDI dataset <https://purl.stanford.edu/ng782rw8378>`_.
+#
+
+afd.organize_stanford_data(clear_previous_afq=True)
+
+##########################################################################
+# Set tractography parameters (optional)
+# ---------------------
+# We make this tracking_params which we will pass to the AFQ object
+# which specifies that we want 100,000 seeds randomly distributed
+# in the ROIs of every bundle.
+#
+# We only do this to make this example faster and consume less space.
+
+tracking_params = dict(seed_mask=RoiMask(),
+                       n_seeds=10000,
+                       random_seeds=True,
+                       rng_seed=42)
+
+##########################################################################
+# Initialize an AFQ object:
+# -------------------------
+#
+# We specify bundle_info as the default bundles list (api.BUNDLES) plus the
+# callosal bundle list. This tells the AFQ object to use bundles from both
+# the standard and callosal templates.
+
+myafq = api.AFQ(bids_path=op.join(afd.afq_home,
+                                  'stanford_hardi'),
+                dmriprep='vistasoft',
+                bundle_info=api.BUNDLES + api.CALLOSUM_BUNDLES,
+                tracking_params=tracking_params)
+
+##########################################################################
+# Visualizing bundles and tract profiles:
+# ---------------------------------------
+# This would run the script and visualize the bundles using the plotly
+# interactive visualization, which should automatically open in a
+# new browser window.
+bundle_html = myafq.viz_bundles(export=True, n_points=50)
+plotly.io.show(bundle_html[0])

examples/plot_afq_reco80.py

@@ -9,8 +9,6 @@
 """
 import os.path as op
 
-import matplotlib.pyplot as plt
-import nibabel as nib
 import plotly
 
 from AFQ import api
@@ -23,7 +21,20 @@
 # Retrieves `Stanford HARDI dataset <https://purl.stanford.edu/ng782rw8378>`_.
 #
 
-afd.organize_stanford_data()
+afd.organize_stanford_data(clear_previous_afq=True)
+
+##########################################################################
+# Set tractography parameters (optional)
+# ---------------------
+# We make this tracking_params which we will pass to the AFQ object
+# which specifies that we want 50,000 seeds randomly distributed
+# in the white matter.
+#
+# We only do this to make this example faster and consume less space.
+
+tracking_params = dict(n_seeds=50000,
+                       random_seeds=True,
+                       rng_seed=42)
 
 ##########################################################################
 # Initialize an AFQ object:
@@ -36,7 +47,8 @@
 myafq = api.AFQ(bids_path=op.join(afd.afq_home,
                                   'stanford_hardi'),
                 dmriprep='vistasoft',
-                segmentation_params={"seg_algo": "reco80"})
+                segmentation_params={"seg_algo": "reco80"},
+                tracking_params=tracking_params)
 
 ##########################################################################
 # Visualizing bundles and tract profiles:

examples/plot_bids_layout.py

@@ -47,7 +47,7 @@
 # this dataset and organized it within the `~/AFQ_data` folder in the BIDS
 # format.
 
-afd.organize_stanford_data()
+afd.organize_stanford_data(clear_previous_afq=True)
 
 ##########################################################################
 # After doing that, we should have a folder that looks like this: