[ENH] Add DSI Studio AutoTrack recon workflow (#576)

.circleci/AutoTrackTest.sh

@@ -0,0 +1,43 @@
+#!/bin/bash
+
+cat << DOC
+
+Reconstruction workflow tests
+=============================
+
+All supported reconstruction workflows get tested
+
+This tests the following features:
+
+Inputs:
+-------
+
+ - qsiprep single shell results (data/DSDTI_fmap)
+ - qsiprep multi shell results (data/DSDTI_fmap)
+
+DOC
+set +e
+source ./get_data.sh
+TESTDIR=${PWD}
+get_config_data ${TESTDIR}
+get_bids_data ${TESTDIR} multishell_output
+CFG=${TESTDIR}/data/nipype.cfg
+EDDY_CFG=${TESTDIR}/data/eddy_config.json
+export FS_LICENSE=${TESTDIR}/data/license.txt
+
+# Test MRtrix3 multishell msmt with ACT
+TESTNAME=autotrack
+setup_dir ${TESTDIR}/${TESTNAME}
+TEMPDIR=${TESTDIR}/${TESTNAME}/work
+OUTPUT_DIR=${TESTDIR}/${TESTNAME}/derivatives
+BIDS_INPUT_DIR=${TESTDIR}/data/multishell_output/qsiprep
+QSIPREP_CMD=$(run_qsiprep_cmd ${BIDS_INPUT_DIR} ${OUTPUT_DIR})
+
+${QSIPREP_CMD} \
+	 -w ${TEMPDIR} \
+	 --recon-input ${BIDS_INPUT_DIR} \
+	 --sloppy \
+     --stop-on-first-crash \
+	 --recon-spec dsi_studio_autotrack \
+	 --recon-only \
+	 -vv

.circleci/config.yml

@@ -197,6 +197,18 @@ jobs:
             cd .circleci
             bash Mrtrix3ReconTests.sh
 
+  Recon_AutoTrack:
+    <<: *dockersetup
+    steps:
+      - checkout
+      - run: *runinstall
+      - run:
+          name: Test the CSD recon workflows
+          no_output_timeout: 1h
+          command: |
+            cd .circleci
+            bash AutoTrackTest.sh
+
   Recon_DIPY:
     <<: *dockersetup
     steps:
@@ -535,6 +547,13 @@ workflows:
             tags:
               only: /.*/
 
+      - Recon_AutoTrack:
+          requires:
+            - build
+          filters:
+            tags:
+              only: /.*/
+
       - Recon_DIPY:
           requires:
             - build
@@ -577,6 +596,7 @@ workflows:
             - MultiT1w
             - Recon_3Tissue
             - Recon_MRtrix3
+            - Recon_AutoTrack
             - Recon_DIPY
             - Recon_AMICO
             - Recon_PYAFQ

qsiprep/data/pipelines/dsi_studio_autotrack.json

@@ -0,0 +1,28 @@
+{
+  "name": "dsistudio_pipeline",
+  "space": "T1w",
+  "atlases": [],
+  "nodes": [
+    {
+      "name": "dsistudio_gqi",
+      "software": "DSI Studio",
+      "action": "reconstruction",
+      "input": "qsiprep",
+      "output_suffix": "gqi",
+      "parameters": {"method": "gqi"}
+    },
+    {
+      "name": "autotrackgqi",
+      "software": "DSI Studio",
+      "action": "autotrack",
+      "input": "dsistudio_gqi",
+      "output_suffix": "AutoTrackGQI",
+      "parameters": {
+        "track_id": "Fasciculus,Cingulum,Aslant,Corticos,Thalamic_R,Reticular,Optic,Fornix,Corpus",
+        "tolerance": "22,26,30",
+        "track_voxel_ratio": 2.0,
+        "yield_rate": 0.000001
+      }
+    }
+  ]
+}

qsiprep/interfaces/bids.py

@@ -276,6 +276,7 @@ class DerivativesDataSinkInputSpec(BaseInterfaceInputSpec):
     source_file = File(mandatory=True, desc='the original file or name of merged files')
     space = traits.Str('', usedefault=True, desc='Label for space field')
     desc = traits.Str('', usedefault=True, desc='Label for description field')
+    bundle = traits.Str('', usedefault=True, desc='Label for bundle field')
     suffix = traits.Str('', usedefault=True, desc='suffix appended to source_file')
     keep_dtype = traits.Bool(False, usedefault=True, desc='keep datatype suffix')
     extra_values = traits.List(traits.Str)
@@ -366,10 +367,12 @@ def _run_interface(self, runtime):
         os.makedirs(out_path, exist_ok=True)
         base_fname = op.join(out_path, src_fname)
 
-        formatstr = '{bname}{space}{desc}{suffix}{dtype}{ext}'
+        formatstr = '{bname}{space}{desc}{bundle}{suffix}{dtype}{ext}'
 
         space = '_space-{}'.format(self.inputs.space) if self.inputs.space else ''
         desc = '_desc-{}'.format(self.inputs.desc) if self.inputs.desc else ''
+        bundle = '_bundle-{}'.format(
+            self.inputs.bundle.replace("_", "")) if self.inputs.bundle else ''
         suffix = '_{}'.format(self.inputs.suffix) if self.inputs.suffix else ''
         dtype = '' if not self.inputs.keep_dtype else ('_%s' % dtype)
 
@@ -380,6 +383,7 @@ def _run_interface(self, runtime):
                 bname=base_fname,
                 space=space,
                 desc=desc,
+                bundle=bundle,
                 suffix=suffix,
                 dtype=dtype,
                 ext='')
@@ -389,7 +393,7 @@ def _run_interface(self, runtime):
             return runtime
 
         if len(self.inputs.in_file) > 1 and not isdefined(self.inputs.extra_values):
-            formatstr = '{bname}{space}{desc}{suffix}{i:04d}{dtype}{ext}'
+            formatstr = '{bname}{space}{desc}{bundle}{suffix}{i:04d}{dtype}{ext}'
 
         # Otherwise it's file(s)
         self._results['compression'] = []
@@ -398,6 +402,7 @@ def _run_interface(self, runtime):
                 bname=base_fname,
                 space=space,
                 desc=desc,
+                bundle=bundle,
                 suffix=suffix,
                 i=i,
                 dtype=dtype,

qsiprep/interfaces/converters.py

@@ -151,6 +151,50 @@ def _run_interface(self, runtime):
         return runtime
 
 
+class _DSIStudioTrkToTckInputSpec(BaseInterfaceInputSpec):
+    trk_file = File(exists=True, mandatory=True)
+    reference_nifti = File(exists=True, mandatory=True)
+
+
+class _DSIStudioTrkToTckOutputSpec(TraitedSpec):
+    tck_file = File()
+
+
+class DSIStudioTrkToTck(SimpleInterface):
+    input_spec = _DSIStudioTrkToTckInputSpec
+    output_spec = _DSIStudioTrkToTckOutputSpec
+
+    def _run_interface(self, runtime):
+
+        if self.inputs.trk_file.endswith(".gz"):
+            with gzip.open(self.inputs.trk_file, "r") as trkf:
+                dsi_trk = nb.streamlines.load(trkf)
+        else:
+            dsi_trk = nb.streamlines.load(self.inputs.trk_file)
+
+        #load preprocessed dwi image
+        dwi_img = nb.load(self.inputs.reference_nifti)
+
+        #convert to voxel coordinates
+        pts = dsi_trk.streamlines._data
+        zooms = np.abs(np.diag(dsi_trk.header['voxel_to_rasmm'])[:3])
+        voxel_coords = pts / zooms
+        voxel_coords[:, 0] = dwi_img.shape[0] - voxel_coords[:, 0]
+        voxel_coords[:, 1] = dwi_img.shape[1] - voxel_coords[:, 1]
+
+        #create new tck
+        new_data = nb.affines.apply_affine(dwi_img.affine, voxel_coords)
+        dsi_trk.tractogram.streamlines._data = new_data
+        tck = nb.streamlines.TckFile(dsi_trk.tractogram)
+        tck_file = fname_presuffix(self.inputs.trk_file.rstrip(".gz"),
+                                   newpath=runtime.cwd,
+                                   use_ext=False,
+                                   suffix=".tck")
+        tck.save(tck_file)
+        self._results['tck_file'] = tck_file
+        return runtime
+
+
 def get_dsi_studio_ODF_geometry(odf_key):
     mat_path = pkgr('qsiprep', 'data/odfs.mat')
     m = loadmat(mat_path)