Merge pull request #417 from dPys/development

Development

Author: dPys

docs/usage.rst

@@ -119,12 +119,12 @@ Command-Line Arguments
 **********************
 
 .. argparse::
-    :module: pynets.cli.pynets_bids
+    :module: pynets.cli.pynets_run
     :func: get_bids_parser
     :prog: pynets
 
 .. argparse::
-    :module: pynets.cli.pynets_run
+    :module: pynets.cli.pynets_bids
     :func: get_parser
     :prog: pynets
 

pynets/cli/pynets_bids.py

@@ -2,6 +2,7 @@
 PyNets BIDS cli
 """
 import bids
+
 from pynets.core.utils import as_list, merge_dicts
 
 
@@ -440,7 +441,6 @@ def main():
     import yaml
     import itertools
     from types import SimpleNamespace
-    from pathlib import Path
     import pkg_resources
     from pynets.core.utils import flatten
     from pynets.cli.pynets_run import build_workflow

pynets/cli/pynets_run.py

@@ -400,6 +400,18 @@ def get_parser():
              "is 20. If you wish to iterate the pipeline across multiple "
              "minimums, separate the list by space (e.g. 10 30 50).\n",
     )
+    parser.add_argument(
+        "-em",
+        metavar="Error margin",
+        default=12,
+        nargs="+",
+        help="(Hyperparameter): Distance (in the units of the streamlines, "
+             "usually mm). If any coordinate in the streamline is within this "
+             "distance from the center of any voxel in the ROI, the filtering "
+             "criterion is set to True for this streamline, otherwise False. "
+             "Defaults to the distance between the center of each voxel and "
+             "the corner of the voxel.\n",
+    )
     parser.add_argument(
         "-dg",
         metavar="Direction getter",
@@ -973,6 +985,20 @@ def build_workflow(args, retval):
             min_length_list = None
     else:
         min_length_list = None
+    error_margin = args.em
+    if error_margin:
+        if (isinstance(error_margin, list)) and (len(error_margin) > 1):
+            error_margin_list = error_margin
+            error_margin = None
+        elif error_margin == ["None"]:
+            error_margin_list = None
+        elif isinstance(error_margin, list):
+            error_margin = error_margin[0]
+            error_margin_list = None
+        else:
+            error_margin_list = None
+    else:
+        error_margin_list = None
     directget = args.dg
     if directget:
         if (isinstance(directget, list)) and (len(directget) > 1):
@@ -1008,19 +1034,21 @@ def build_workflow(args, retval):
     ) as stream:
         try:
             hardcoded_params = yaml.load(stream)
-            maxcrossing = hardcoded_params["maxcrossing"][0]
+            maxcrossing = hardcoded_params['tracking']["maxcrossing"][0]
             local_corr = hardcoded_params["clustering_local_conn"][0]
-            track_type = hardcoded_params["tracking_method"][0]
-            tiss_class = hardcoded_params["tissue_classifier"][0]
-            target_samples = hardcoded_params["tracking_samples"][0]
+            track_type = hardcoded_params['tracking']["tracking_method"][0]
+            tiss_class = hardcoded_params['tracking']["tissue_classifier"][0]
+            target_samples = hardcoded_params['tracking']["tracking_samples"][0]
             use_parcel_naming = hardcoded_params["parcel_naming"][0]
-            step_list = hardcoded_params["step_list"]
-            curv_thr_list = hardcoded_params["curv_thr_list"]
+            step_list = hardcoded_params['tracking']["step_list"]
+            curv_thr_list = hardcoded_params['tracking']["curv_thr_list"]
             nilearn_parc_atlases = hardcoded_params["nilearn_parc_atlases"]
             nilearn_coord_atlases = hardcoded_params["nilearn_coord_atlases"]
             nilearn_prob_atlases = hardcoded_params["nilearn_prob_atlases"]
             local_atlases = hardcoded_params["local_atlases"]
             template_name = hardcoded_params['template'][0]
+            roi_neighborhood_tol = \
+                hardcoded_params['tracking']["roi_neighborhood_tol"][0]
 
             if track_type == "particle":
                 tiss_class = "cmc"
@@ -1800,6 +1828,24 @@ def build_workflow(args, retval):
         else:
             print(f"{Fore.GREEN}Iterating minimum streamline lengths:")
             print(f"{Fore.BLUE}{', '.join(min_length_list)}")
+        if error_margin:
+            if float(roi_neighborhood_tol) <= float(error_margin):
+                raise ValueError(
+                    'roi_neighborhood_tol preset cannot be less than '
+                    'the value of the structural connectome error_margin'
+                    ' parameter.')
+            print(f"{Fore.GREEN}Using {Fore.BLUE}{error_margin}"
+                  f"mm{Fore.GREEN} error margin...")
+        else:
+            for em in error_margin_list:
+                if float(roi_neighborhood_tol) <= float(em):
+                    raise ValueError(
+                        'roi_neighborhood_tol preset cannot be less than '
+                        'the value of the structural connectome error_margin'
+                        ' parameter.')
+            print(f"{Fore.GREEN}Iterating minimum streamline lengths:")
+            print(f"{Fore.BLUE}{', '.join(error_margin_list)}")
+
         if target_samples:
             print(f"{Fore.GREEN}Using {Fore.BLUE}{target_samples} "
                   f"{Fore.GREEN}streamline samples...")
@@ -2217,6 +2263,7 @@ def init_wf_single_subject(
         track_type,
         min_length,
         maxcrossing,
+        error_margin,
         directget,
         tiss_class,
         runtime_dict,
@@ -2231,6 +2278,7 @@ def init_wf_single_subject(
         waymask,
         local_corr,
         min_length_list,
+        error_margin_list,
         extract_strategy,
         extract_strategy_list,
         outdir,
@@ -2384,6 +2432,7 @@ def init_wf_single_subject(
                 track_type,
                 min_length,
                 maxcrossing,
+                error_margin,
                 directget,
                 tiss_class,
                 runtime_dict,
@@ -2398,6 +2447,7 @@ def init_wf_single_subject(
                 waymask,
                 local_corr,
                 min_length_list,
+                error_margin_list,
                 extract_strategy,
                 extract_strategy_list,
                 outdir,
@@ -2650,6 +2700,7 @@ def wf_multi_subject(
         track_type,
         min_length,
         maxcrossing,
+        error_margin,
         directget,
         tiss_class,
         runtime_dict,
@@ -2664,6 +2715,7 @@ def wf_multi_subject(
         waymask,
         local_corr,
         min_length_list,
+        error_margin_list,
         extract_strategy,
         extract_strategy_list,
         outdir,
@@ -2776,6 +2828,7 @@ def wf_multi_subject(
                     track_type=track_type,
                     min_length=min_length,
                     maxcrossing=maxcrossing,
+                    error_margin=error_margin,
                     directget=directget,
                     tiss_class=tiss_class,
                     runtime_dict=runtime_dict,
@@ -2790,6 +2843,7 @@ def wf_multi_subject(
                     waymask=waymask,
                     local_corr=local_corr,
                     min_length_list=min_length_list,
+                    error_margin_list=error_margin_list,
                     extract_strategy=extract_strategy,
                     extract_strategy_list=extract_strategy_list,
                     outdir=subj_dir,
@@ -2893,6 +2947,7 @@ def wf_multi_subject(
                     track_type=track_type,
                     min_length=min_length,
                     maxcrossing=maxcrossing,
+                    error_margin=error_margin,
                     directget=directget,
                     tiss_class=tiss_class,
                     runtime_dict=runtime_dict,
@@ -2907,6 +2962,7 @@ def wf_multi_subject(
                     waymask=waymask,
                     local_corr=local_corr,
                     min_length_list=min_length_list,
+                    error_margin_list=error_margin_list,
                     extract_strategy=extract_strategy,
                     extract_strategy_list=extract_strategy_list,
                     outdir=subj_dir,
@@ -3008,6 +3064,7 @@ def wf_multi_subject(
             track_type,
             min_length,
             maxcrossing,
+            error_margin,
             directget,
             tiss_class,
             runtime_dict,
@@ -3022,6 +3079,7 @@ def wf_multi_subject(
             waymask,
             local_corr,
             min_length_list,
+            error_margin_list,
             extract_strategy,
             extract_strategy_list,
             outdir,
@@ -3195,6 +3253,7 @@ def wf_multi_subject(
             track_type,
             min_length,
             maxcrossing,
+            error_margin,
             directget,
             tiss_class,
             runtime_dict,
@@ -3209,6 +3268,7 @@ def wf_multi_subject(
             waymask,
             local_corr,
             min_length_list,
+            error_margin_list,
             extract_strategy,
             extract_strategy_list,
             subj_dir,

pynets/core/interfaces.py

@@ -952,7 +952,7 @@ class _PlotStructInputSpec(BaseInterfaceInputSpec):
     track_type = traits.Any(mandatory=True)
     directget = traits.Any(mandatory=True)
     min_length = traits.Any(mandatory=True)
-
+    error_margin = traits.Any(mandatory=True)
 
 class _PlotStructOutputSpec(BaseInterfaceInputSpec):
     """Output interface wrapper for PlotStruct"""
@@ -998,6 +998,7 @@ def _run_interface(self, runtime):
                 self.inputs.track_type,
                 self.inputs.directget,
                 self.inputs.min_length,
+                self.inputs.error_margin
             )
 
         self._results["out"] = "None"
@@ -2525,6 +2526,7 @@ class _TrackingInputSpec(BaseInterfaceInputSpec):
     step_list = traits.List(mandatory=True)
     track_type = traits.Str(mandatory=True)
     min_length = traits.Any(mandatory=True)
+    error_margin = traits.Any(mandatory=True)
     maxcrossing = traits.Any(mandatory=True)
     directget = traits.Str(mandatory=True)
     conn_model = traits.Str(mandatory=True)
@@ -2549,8 +2551,6 @@ class _TrackingInputSpec(BaseInterfaceInputSpec):
     fa_path = File(exists=True, mandatory=True)
     waymask = traits.Any(mandatory=False)
     t1w2dwi = File(exists=True, mandatory=True)
-    roi_neighborhood_tol = traits.Any(10, mandatory=True, usedefault=True)
-    sphere = traits.Str('repulsion724', mandatory=True, usedefault=True)
 
 
 class _TrackingOutputSpec(TraitedSpec):
@@ -2583,9 +2583,8 @@ class _TrackingOutputSpec(TraitedSpec):
     dm_path = File(exists=True, mandatory=True)
     directget = traits.Str(mandatory=True)
     labels_im_file = File(exists=True, mandatory=True)
-    roi_neighborhood_tol = traits.Any()
     min_length = traits.Any()
-
+    error_margin = traits.Any()
 
 class Tracking(SimpleInterface):
     """Interface wrapper for Tracking"""
@@ -2617,7 +2616,10 @@ def _run_interface(self, runtime):
             pkg_resources.resource_filename("pynets", "runconfig.yaml"), "r"
         ) as stream:
             hardcoded_params = yaml.load(stream)
-            use_life = hardcoded_params["use_life"][0]
+            use_life = hardcoded_params['tracking']["use_life"][0]
+            roi_neighborhood_tol = hardcoded_params['tracking']["roi_neighborhood_tol"][0]
+            sphere = hardcoded_params['tracking']["sphere"][0]
+
         stream.close()
 
         dir_path = utils.do_dir_path(
@@ -2682,12 +2684,12 @@ def _run_interface(self, runtime):
             "_",
             "%s"
             % (
-                "%s%s" % (self.inputs.node_size, "mm_")
+                "%s%s" % (self.inputs.node_size, "mm")
                 if (
                     (self.inputs.node_size != "parc")
                     and (self.inputs.node_size is not None)
                 )
-                else "parc_"
+                else "parc"
             ),
             ".npy",
         )
@@ -2829,16 +2831,17 @@ def _run_interface(self, runtime):
 
         # Iteratively build a list of streamlines for each ROI while tracking
         print(
-            f"{Fore.GREEN}Target number of samples: {Fore.BLUE} "
+            f"{Fore.GREEN}Target number of cumulative streamlines: "
+            f"{Fore.BLUE} "
             f"{self.inputs.target_samples}"
         )
         print(Style.RESET_ALL)
         print(
-            f"{Fore.GREEN}Using curvature threshold(s): {Fore.BLUE} "
+            f"{Fore.GREEN}Curvature threshold(s): {Fore.BLUE} "
             f"{self.inputs.curv_thr_list}"
         )
         print(Style.RESET_ALL)
-        print(f"{Fore.GREEN}Using step size(s): {Fore.BLUE} "
+        print(f"{Fore.GREEN}Step size(s): {Fore.BLUE} "
               f"{self.inputs.step_list}")
         print(Style.RESET_ALL)
         print(f"{Fore.GREEN}Tracking type: {Fore.BLUE} "
@@ -2865,13 +2868,13 @@ def _run_interface(self, runtime):
             atlas_data_wm_gm_int,
             parcels,
             model,
-            get_sphere(self.inputs.sphere),
+            get_sphere(sphere),
             self.inputs.directget,
             self.inputs.curv_thr_list,
             self.inputs.step_list,
             self.inputs.track_type,
             self.inputs.maxcrossing,
-            int(self.inputs.roi_neighborhood_tol),
+            int(roi_neighborhood_tol),
             self.inputs.min_length,
             waymask_data,
             B0_mask_data,
@@ -2880,7 +2883,7 @@ def _run_interface(self, runtime):
             self.inputs.tiss_class, B0_mask_tmp_path
         )
 
-        del model, parcels, atlas_data_wm_gm_int, waymask_data
+        del model, parcels, atlas_data_wm_gm_int
         gc.collect()
 
         # Save streamlines to trk
@@ -2927,10 +2930,16 @@ def _run_interface(self, runtime):
             dwi_img = nib.load(dwi_file_tmp_path, mmap=True)
             dwi_data = dwi_img.get_fdata().astype('float32')
             orig_count = len(streamlines)
+
+            if self.inputs.waymask:
+                mask_data = waymask_data
+            else:
+                mask_data = nib.load(wm_in_dwi_tmp_path
+                                     ).get_fdata().astype('float32')
             try:
                 streamlines = evaluate_streamline_plausibility(
-                    dwi_data, gtab, B0_mask_data, streamlines,
-                    sphere=self.inputs.sphere)
+                    dwi_data, gtab, mask_data, streamlines,
+                    sphere=sphere)
             except BaseException:
                 print(f"Linear Fascicle Evaluation failed. Visually checking "
                       f"streamlines output {namer_dir}/{op.basename(streams)}"
@@ -2940,7 +2949,7 @@ def _run_interface(self, runtime):
                                  'the tractogram!')
             del dwi_data
 
-        del B0_mask_data
+        del B0_mask_data, waymask_data
 
         stf = StatefulTractogram(
             streamlines,
@@ -3011,9 +3020,8 @@ def _run_interface(self, runtime):
         self._results["dm_path"] = dm_path
         self._results["directget"] = self.inputs.directget
         self._results["labels_im_file"] = labels_im_file_tmp_path
-        self._results["roi_neighborhood_tol"] = \
-            self.inputs.roi_neighborhood_tol
         self._results["min_length"] = self.inputs.min_length
+        self._results["error_margin"] = self.inputs.error_margin
 
         tmp_files = [B0_mask_tmp_path, gtab_file_tmp_path,
                      labels_im_file_tmp_path_wm_gm_int, dwi_file_tmp_path]