[ENH] Update DSI Studio to the latest commit (#573)

.circleci/config.yml

@@ -2,7 +2,7 @@ version: 2.1
 
 .dockersetup: &dockersetup
   docker:
-    - image: pennbbl/qsiprep_build:23.5.5
+    - image: pennbbl/qsiprep_build:23.5.6
   working_directory: /src/qsiprep
 
 runinstall: &runinstall

Dockerfile

@@ -1,4 +1,4 @@
-FROM pennbbl/qsiprep_build:23.5.5
+FROM pennbbl/qsiprep_build:23.5.6
 
 # WORKDIR /root/
 # # Installing qsiprep

qsiprep/data/pipelines/dsi_studio_gqi.json

@@ -25,7 +25,6 @@
       "input": "dsistudio_gqi",
       "parameters": {
         "turning_angle": 35,
-        "method": 0,
         "smoothing": 0.0,
         "step_size": 1.0,
         "min_length": 30,

qsiprep/interfaces/dsi_studio.py

@@ -1,5 +1,6 @@
 #!python
 import os
+from pathlib import Path
 import os.path as op
 import logging
 from glob import glob
@@ -34,6 +35,10 @@ class DSIStudioCreateSrcInputSpec(DSIStudioCommandLineInputSpec):
         desc="Directory with DICOM data from only the dwi",
         exists=True,
         argstr="--source=%s")
+    bvec_convention = traits.Enum(
+        ("DIPY", "FSL"),
+        usedefault=True,
+        desc="Convention used for bvecs. FSL assumes LAS+ no matter image orientation")
     input_bvals_file = File(
         desc="Text file containing b values", exists=True, argstr="--bval=%s")
     input_bvecs_file = File(
@@ -72,6 +77,32 @@ class DSIStudioCreateSrc(CommandLine):
     output_spec = DSIStudioCreateSrcOutputSpec
     _cmd = "dsi_studio --action=src "
 
+    def _pre_run_hook(self, runtime):
+        """As of QSIPrep > 0.17 DSI Studio changed from DIPY bvecs to FSL bvecs."""
+
+        # b_table files and dicom directories are ok
+        if isdefined(self.inputs.input_b_table_file) or \
+                isdefined(self.inputs.input_dicom_dir):
+            return runtime
+
+        if not (isdefined(self.inputs.input_bvals_file) and \
+                isdefined(self.inputs.input_bvecs_file)):
+            raise Exception("without a b_table or dicom directory, both bvals and bvecs "
+                            "must be specified")
+
+        # If the bvecs are in DIPY format, convert them to a b_table.txt
+        if self.inputs.bvec_convention == "DIPY":
+            btable_file = self._gen_filename('output_src').replace(".src.gz", ".b_table.txt")
+            btable_from_bvals_bvecs(self.inputs.input_bvals_file,
+                                    self.inputs.input_bvecs_file,
+                                    btable_file)
+            self.inputs.input_b_table_file = btable_file
+            self.inputs.input_bvals_file = traits.Undefined
+            self.inputs.input_bvecs_file = traits.Undefined
+            LOGGER.info("Converted DIPY LPS+ bval/bvec to DSI Studio b_table")
+        return runtime
+
+
     def _gen_filename(self, name):
         if not name == 'output_src':
             return None
@@ -155,6 +186,19 @@ class DSIStudioReconstructionInputSpec(DSIStudioCommandLineInputSpec):
         position=-1,
         argstr="#%s")
     thread_count = traits.Int(1, usedefault=True, argstr="--thread_count=%d")
+
+    dti_no_high_b = traits.Bool(
+        True,
+        usedefault=True,
+        argstr="--dti_no_high_b=%d",
+        desc="specify whether the construction of DTI should ignore high b-value (b>1500)")
+    r2_weighted = traits.Bool(
+        False,
+        usedefault=True,
+        argstr="--r2_weighted=%d",
+        desc="specify whether GQI and QSDR uses r2-weighted to calculate SDF")
+
+    # Outputs
     output_odf = traits.Bool(
         True,
         usedefault=True,
@@ -163,40 +207,31 @@ class DSIStudioReconstructionInputSpec(DSIStudioCommandLineInputSpec):
     odf_order = traits.Enum((8, 4, 5, 6, 10, 12, 16, 20),
                             usedefault=True,
                             desc="ODF tesselation order")
+    # Which scalars to include
+    other_output = traits.Str(
+        "all",
+        argstr="--other_output=%s",
+        desc="additional diffusion metrics to calculate",
+        usedefault=True
+    )
+    align_acpc = traits.Bool(
+        False,
+        usedefault=True,
+        argstr="--align_acpc=%d",
+        desc="rotate image volume to align ap-pc"
+    )
     check_btable = traits.Enum(
         (0, 1),
         usedefault=True,
         argstr="--check_btable=%d",
         desc="Check if btable matches nifti orientation (not foolproof)")
+
     num_fibers = traits.Int(
         3,
         usedefault=True,
         argstr="--num_fiber=%d",
         desc="number of fiber populations estimated at each voxel")
 
-    # Decomposition traits
-    decomposition = traits.Bool(
-        False, argstr="--decomposition=1", desc="Apply ODF Decomposition")
-    decomp_fraction = traits.Float(
-        0.05,
-        desc="Decomposition Fraction",
-        argstr="--param3=%.2f",
-        requires=["decomposition"])
-    decomp_m_value = traits.Int(
-        10,
-        desc="Decomposition m value",
-        argstr="--param4=%d",
-        requires=["decomposition"])
-
-    # Deconvolution traits
-    deconvolution = traits.Bool(
-        False, argstr="--deconvolution=1", desc="Apply ODF Deconvolution")
-    deconv_regularization = traits.Float(
-        0.5,
-        requires=["deconvolution"],
-        desc="Deconvolution regularization parameter",
-        argstr="--param2=%.2f")
-
 
 class DSIStudioReconstructionOutputSpec(TraitedSpec):
     output_fib = File(desc="Output File", exists=True)
@@ -262,17 +297,27 @@ class DSIStudioExportInputSpec(DSIStudioCommandLineInputSpec):
 
 
 class DSIStudioExportOutputSpec(DSIStudioCommandLineInputSpec):
-    gfa_file = File(desc="Exported files")
-    fa0_file = File(desc="Exported files")
-    fa1_file = File(desc="Exported files")
-    fa2_file = File(desc="Exported files")
-    fa3_file = File(desc="Exported files")
-    fa4_file = File(desc="Exported files")
-    iso_file = File(desc="Exported files")
-    dti_fa_file = File(desc="Exported files")
-    md_file = File(desc="Exported files")
-    rd_file = File(desc="Exported files")
-    ad_file = File(desc="Exported files")
+    qa_file = File(desc="Exported scalar nifti")
+    color_file = File(desc="Exported scalar nifti")
+    dti_fa_file = File(desc="Exported scalar nifti")
+    txx_file = File(desc="Exported scalar nifti")
+    txy_file = File(desc="Exported scalar nifti")
+    txz_file = File(desc="Exported scalar nifti")
+    tyy_file = File(desc="Exported scalar nifti")
+    tyz_file = File(desc="Exported scalar nifti")
+    tzz_file = File(desc="Exported scalar nifti")
+    rd1_file = File(desc="Exported scalar nifti")
+    rd2_file = File(desc="Exported scalar nifti")
+    ha_file = File(desc="Exported scalar nifti")
+    md_file = File(desc="Exported scalar nifti")
+    ad_file = File(desc="Exported scalar nifti")
+    rd_file = File(desc="Exported scalar nifti")
+    gfa_file = File(desc="Exported scalar nifti")
+    iso_file = File(desc="Exported scalar nifti")
+    rdi_file = File(desc="Exported scalar nifti")
+    nrdi02L_file = File(desc="Exported scalar nifti")
+    nrdi04L_file = File(desc="Exported scalar nifti")
+    nrdi06L_file = File(desc="Exported scalar nifti")
     image0_file = File(desc="Exported files")
 
 
@@ -284,22 +329,29 @@ class DSIStudioExport(CommandLine):
     def _list_outputs(self):
         outputs = self.output_spec().get()
         to_expect = self.inputs.to_export.split(",")
+        cwd = Path()
+        results = list(cwd.glob("*.nii.gz"))
         for expected in to_expect:
-            matches = glob("*" + expected + "*.nii.gz")
+            matches = [fname.absolute().name for fname in results if
+                       fname.name.endswith("." + expected + ".nii.gz")]
             if len(matches) == 1:
-                outputs[expected + "_file"] = os.path.abspath(matches[0])
+                outputs[expected + "_file"] = matches[0]
+            elif len(matches) == 0:
+                raise Exception("No exported scalar found matching " + expected)
+            else:
+                raise Exception("Found too mand scalar files matching " + expected)
+
         return outputs
 
 
 class DSIStudioConnectivityMatrixInputSpec(DSIStudioCommandLineInputSpec):
-    trk_file = File(exists=True, argstr="--tract=%s")
+    trk_file = File(exists=True, argstr="--tract=%s", copyfile=False)
     atlas_config = traits.Dict(desc='atlas configs for atlases to run connectivity for')
     atlas_name = traits.Str()
     input_fib = File(
         exists=True, argstr="--source=%s", mandatory=True, copyfile=False)
     fiber_count = traits.Int(xor=["seed_count"], argstr="--fiber_count=%d")
     seed_count = traits.Int(xor=["fiber_count"], argstr="--seed_count=%d")
-    method = traits.Enum((0, 4), argstr="--method=%d", usedefault=True)
     seed_plan = traits.Enum((0, 1), argstr="--seed_plan=%d")
     initial_dir = traits.Enum((0, 1, 2), argstr="--initial_dir=%d")
     interpolation = traits.Enum((0, 1, 2), argstr="--interpolation=%d")
@@ -640,7 +692,7 @@ def _run_interface(self, runtime):
 
         # No matter what, still use the correct affine
         nb.Nifti1Image(
-            reoriented_img.get_fdata()[::-1, ::-1, :],
+            reoriented_img.get_fdata(),
             correct_img.affine).to_filename(new_file)
         self._results['out_file'] = new_file
 
@@ -671,6 +723,32 @@ def _run_interface(self, runtime):
         return runtime
 
 
+class _DSIStudioBTableInputSpec(BaseInterfaceInputSpec):
+    bval_file = File(exists=True, mandatory=True)
+    bvec_file = File(exists=True, mandatory=True)
+    bvec_convention = traits.Enum(
+        ("DIPY", "FSL"),
+        usedefault=True,
+        desc="Convention used for bvecs. FSL assumes LAS+ no matter image orientation")
+
+
+class _DSIStudioBTableOutputSpec(TraitedSpec):
+    btable_file = File(exists=True)
+
+
+class DSIStudioBTable(SimpleInterface):
+    input_spec = _DSIStudioBTableInputSpec
+    output_spec = _DSIStudioBTableOutputSpec
+
+    def _run_interface(self, runtime):
+        if self.inputs.bvec_convention != "DIPY":
+            raise NotImplementedError("Only DIPY Bvecs supported for now")
+        btab_file = op.join(runtime.cwd, "btable.txt")
+        btable_from_bvals_bvecs(self.inputs.bval_file, self.inputs.bvec_file, btab_file)
+        self._results['btable_file'] = btab_file
+        return runtime
+
+
 def load_src_qc_file(fname, prefix=""):
     with open(fname, "r") as qc_file:
         qc_data = qc_file.readlines()
@@ -709,3 +787,29 @@ def load_fib_qc_file(fname):
         lines = [line.strip().split() for line in fibqc_f]
     return {'coherence_index': [float(lines[0][-1])],
             'incoherence_index': [float(lines[1][-1])]}
+
+
+def btable_from_bvals_bvecs(bval_file, bvec_file, output_file):
+    """Create a b-table from DIPY-style bvals/bvecs.
+
+    Assuming these come from qsiprep they will be in LPS+, which
+    is the same convention as DSI Studio's btable.
+    """
+    bvals = np.loadtxt(bval_file).squeeze()
+    bvecs = np.loadtxt(bvec_file).squeeze()
+    if not 3 in bvecs.shape:
+        raise Exception(
+            "uninterpretable bval/bvec files\n\t{}\n\t{}".format(bval_file, bvec_file))
+    if not bvecs.shape[1] == 3:
+        bvecs = bvecs.T
+
+    if not bvecs.shape[0] == bvals.shape[0]:
+        raise Exception("Bval/Bvec mismatch")
+
+    rows = []
+    for row in map(tuple, np.column_stack([bvals, bvecs])):
+        rows.append("%d %.6f %.6f %.6f" % row)
+
+    # Write the actual file:
+    with open(output_file, "w") as btablef:
+        btablef.write("\n".join(rows) + "\n")

qsiprep/workflows/dwi/derivatives.py

@@ -32,7 +32,7 @@ def init_dwi_derivatives_wf(output_prefix,
     inputnode = pe.Node(
         niu.IdentityInterface(fields=[
             'source_file', 'dwi_t1', 'dwi_mask_t1', 'cnr_map_t1', 'bvals_t1', 'bvecs_t1',
-            'local_bvecs_t1', 't1_b0_ref', 'gradient_table_t1', 'confounds',
+            'local_bvecs_t1', 't1_b0_ref', 'gradient_table_t1', 'btable_t1', 'confounds',
             'hmc_optimization_data', 'series_qc'
         ]),
         name='inputnode')
@@ -130,6 +130,17 @@ def init_dwi_derivatives_wf(output_prefix,
         name='ds_gradient_table_t1',
         run_without_submitting=True,
         mem_gb=DEFAULT_MEMORY_MIN_GB)
+    ds_btable_t1 = pe.Node(
+        DerivativesDataSink(
+            source_file=source_file,
+            base_directory=output_dir,
+            space='T1w',
+            suffix='dwi',
+            extension='.b_table.txt',
+            desc='preproc'),
+        name='ds_btable_t1',
+        run_without_submitting=True,
+        mem_gb=DEFAULT_MEMORY_MIN_GB)
 
     workflow.connect([
         (inputnode, ds_dwi_t1, [('dwi_t1', 'in_file')]),
@@ -138,7 +149,8 @@ def init_dwi_derivatives_wf(output_prefix,
         (inputnode, ds_t1_b0_ref, [('t1_b0_ref', 'in_file')]),
         (inputnode, ds_dwi_mask_t1, [('dwi_mask_t1', 'in_file')]),
         (inputnode, ds_cnr_map_t1, [('cnr_map_t1', 'in_file')]),
-        (inputnode, ds_gradient_table_t1, [('gradient_table_t1', 'in_file')])
+        (inputnode, ds_gradient_table_t1, [('gradient_table_t1', 'in_file')]),
+        (inputnode, ds_btable_t1, [('btable_t1', 'in_file')]),
         ])
     # If requested, write local bvecs
     if write_local_bvecs:

qsiprep/workflows/dwi/distortion_group_merge.py

@@ -150,8 +150,12 @@ def init_distortion_group_merge_wf(merging_strategy, inputs_list, hmc_model, rep
     ])
 
     # Calculate QC on the merged raw and processed data
-    raw_qc_wf = init_modelfree_qc_wf(omp_nthreads=omp_nthreads, name='raw_qc_wf')
-    processed_qc_wf = init_modelfree_qc_wf(omp_nthreads=omp_nthreads, name='processed_qc_wf')
+    raw_qc_wf = init_modelfree_qc_wf(omp_nthreads=omp_nthreads,
+                                     bvec_convention="auto",
+                                     name='raw_qc_wf')
+    processed_qc_wf = init_modelfree_qc_wf(omp_nthreads=omp_nthreads,
+                                           bvec_convention="DIPY",  # Always LPS+ when resampled
+                                           name='processed_qc_wf')
     # Combine all the QC measures for a series QC
     series_qc = pe.Node(SeriesQC(output_file_name=output_prefix), name='series_qc')
     ds_series_qc = pe.Node(
@@ -262,6 +266,7 @@ def init_distortion_group_merge_wf(merging_strategy, inputs_list, hmc_model, rep
             ('local_bvecs_t1', 'inputnode.local_bvecs_t1'),
             ('t1_b0_ref', 'inputnode.t1_b0_ref'),
             ('gradient_table_t1', 'inputnode.gradient_table_t1'),
+            ('btable_t1', 'inputnode.btable_t1'),
             ('confounds', 'inputnode.confounds'),
             ('hmc_optimization_data', 'inputnode.hmc_optimization_data')]),
     ])