ENH: Move from DTI to DKI with multishell data

mriqc/data/bootstrap-dwi.yml

@@ -38,11 +38,11 @@ sections:
       the signal variability, and therefore should be interpreted with care.
     subtitle: Shell-wise joint distribution of SNR vs. FA in every voxel
   - bids: {datatype: figures, desc: fa}
-    caption: Reconstructed FA map using Dipy's free-water DTI model.
+    caption: Reconstructed FA map.
     subtitle: Fractional anisotropy (FA) map
-  - bids: {datatype: figures, desc: ad}
-    caption: Reconstructed ADC map using Dipy's free-water DTI model.
-    subtitle: Axial diffusivity (AD) map
+  - bids: {datatype: figures, desc: md}
+    caption: Reconstructed MD map.
+    subtitle: Mean diffusivity (MD) map
 - name: DWI shells
   ordering: bval
   reportlets:

mriqc/interfaces/diffusion.py

@@ -52,7 +52,7 @@
     'CCSegmentation',
     'CorrectSignalDrift',
     'DiffusionQC',
-    'DipyDTI',
+    'DiffusionModel',
     'ExtractOrientations',
     'FilterShells',
     'NumberOfShells',
@@ -697,17 +697,16 @@ def _run_interface(self, runtime):
         return runtime
 
 
-class _DipyDTIInputSpec(_BaseInterfaceInputSpec):
+class _DiffusionModelInputSpec(_BaseInterfaceInputSpec):
     in_file = File(exists=True, mandatory=True, desc='dwi file')
     bvals = traits.List(traits.Float, mandatory=True, desc='bval table')
     bvec_file = File(exists=True, mandatory=True, desc='b-vectors')
-    brainmask = File(exists=True, desc='brain mask file')
-    free_water_model = traits.Bool(False, usedefault=True, desc='use free water model')
-    b_threshold = traits.Float(1100, usedefault=True, desc='use only inner shells of the data')
+    brain_mask = File(exists=True, desc='brain mask file')
     decimals = traits.Int(3, usedefault=True, desc='round output maps for reliability')
+    n_shells = traits.Int(mandatory=True, desc='number of shells')
 
 
-class _DipyDTIOutputSpec(_TraitedSpec):
+class _DiffusionModelOutputSpec(_TraitedSpec):
     out_fa = File(exists=True, desc='output FA file')
     out_fa_nans = File(exists=True, desc='binary mask of NaN values in the "raw" FA map')
     out_fa_degenerate = File(
@@ -718,44 +717,53 @@ class _DipyDTIOutputSpec(_TraitedSpec):
     out_md = File(exists=True, desc='output MD file')
 
 
-class DipyDTI(SimpleInterface):
-    """Split a DWI dataset into ."""
+class DiffusionModel(SimpleInterface):
+    """
+    Fit a :obj:`~dipy.reconst.dki.DiffusionKurtosisModel` on the dataset.
+
+    If ``n_shells`` is set to 1, then a :obj:`~dipy.reconst.dti.TensorModel`
+    is used.
 
-    input_spec = _DipyDTIInputSpec
-    output_spec = _DipyDTIOutputSpec
+    """
+
+    input_spec = _DiffusionModelInputSpec
+    output_spec = _DiffusionModelOutputSpec
 
     def _run_interface(self, runtime):
         from dipy.core.gradients import gradient_table_from_bvals_bvecs
-        from dipy.reconst.dti import TensorModel, color_fa, fractional_anisotropy
-        from dipy.reconst.fwdti import FreeWaterTensorModel
         from nipype.utils.filemanip import fname_presuffix
 
         bvals = np.array(self.inputs.bvals)
-        bval_mask = bvals < self.inputs.b_threshold
 
         gtab = gradient_table_from_bvals_bvecs(
-            bvals=bvals[bval_mask],
-            bvecs=np.loadtxt(self.inputs.bvec_file).T[bval_mask],
+            bvals=bvals,
+            bvecs=np.loadtxt(self.inputs.bvec_file).T,
         )
 
         img = nb.load(self.inputs.in_file)
-        data = img.get_fdata(dtype='float32')[..., bval_mask]
+        data = img.get_fdata(dtype='float32')
 
         brainmask = np.ones_like(data[..., 0], dtype=bool)
 
-        if isdefined(self.inputs.brainmask):
-            brainmask = np.asanyarray(nb.load(self.inputs.brainmask).dataobj) > 0.5
+        if isdefined(self.inputs.brain_mask):
+            brainmask = np.round(
+                nb.load(self.inputs.brain_mask).get_fdata(),
+                3,
+            ) > 0.5
 
-        DTIModel = FreeWaterTensorModel if self.inputs.free_water_model else TensorModel
+        if self.inputs.n_shells == 1:
+            from dipy.reconst.dti import TensorModel as Model
+        else:
+            from dipy.reconst.dki import DiffusionKurtosisModel as Model
 
         # Fit DIT
-        fwdtifit = DTIModel(gtab).fit(
+        fwdtifit = Model(gtab).fit(
             data,
             mask=brainmask,
         )
 
         # Extract the FA
-        fa_data = fractional_anisotropy(fwdtifit.evals)
+        fa_data = fwdtifit.fa
         fa_nan_msk = np.isnan(fa_data)
         fa_data[fa_nan_msk] = 0
 
@@ -823,7 +831,7 @@ def _run_interface(self, runtime):
         fa_degenerate_nii.to_filename(self._results['out_fa_degenerate'])
 
         # Extract the color FA
-        cfa_data = color_fa(fa_data, fwdtifit.evecs)
+        cfa_data = fwdtifit.color_fa
         cfa_nii = nb.Nifti1Image(
             np.clip(cfa_data, a_min=0.0, a_max=1.0),
             img.affine,
@@ -848,15 +856,15 @@ def _run_interface(self, runtime):
             suffix='md',
             newpath=runtime.cwd,
         )
-        ad_data = np.array(fwdtifit.ad, dtype='float32')
-        ad_data[np.isnan(ad_data)] = 0
-        ad_data = np.clip(ad_data, 0, 1)
-        ad_hdr = fa_nii.header.copy()
-        ad_hdr.set_intent('estimate', name='Mean diffusivity (MD)')
+        md_data = np.array(fwdtifit.md, dtype='float32')
+        md_data[np.isnan(md_data)] = 0
+        md_data = np.clip(md_data, 0, 1)
+        md_hdr = fa_nii.header.copy()
+        md_hdr.set_intent('estimate', name='Mean diffusivity (MD)')
         nb.Nifti1Image(
-            ad_data,
+            md_data,
             img.affine,
-            ad_hdr
+            md_hdr
         ).to_filename(self._results['out_md'])
 
         return runtime

mriqc/workflows/diffusion/base.py

@@ -74,9 +74,8 @@ def dmri_qc_workflow(name='dwiMRIQC'):
         PIESNO,
         CCSegmentation,
         CorrectSignalDrift,
-        DipyDTI,
+        DiffusionModel,
         ExtractOrientations,
-        FilterShells,
         NumberOfShells,
         ReadDWIMetadata,
         SpikingVoxelsMask,
@@ -176,8 +175,6 @@ def dmri_qc_workflow(name='dwiMRIQC'):
         iterfield=['in_weights'],
     )
 
-    # Fit DTI model
-    dti_filter = pe.Node(FilterShells(), name='dti_filter')
     dwidenoise = pe.Node(
         DWIDenoise(
             noise='noisemap.nii.gz',
@@ -186,10 +183,8 @@ def dmri_qc_workflow(name='dwiMRIQC'):
         name='dwidenoise',
         nprocs=config.nipype.omp_nthreads,
     )
-    dti = pe.Node(
-        DipyDTI(free_water_model=False),
-        name='dti',
-    )
+    # Fit DTI/DKI model
+    dwimodel = pe.Node(DiffusionModel(), name='dwimodel')
 
     sp_mask = pe.Node(SpikingVoxelsMask(), name='sp_mask')
 
@@ -240,30 +235,28 @@ def dmri_qc_workflow(name='dwiMRIQC'):
         (shells, averages, [('b_masks', 'in_weights')]),
         (averages, hmcwf, [(('out_file', _first), 'inputnode.reference')]),
         (shells, stddev, [('b_masks', 'in_weights')]),
-        (shells, dti_filter, [('out_data', 'bvals')]),
-        (meta, dti_filter, [('out_bvec_file', 'bvec_file')]),
-        (drift, dti_filter, [('out_full_file', 'in_file')]),
-        (dti_filter, dti, [('out_bvals', 'bvals')]),
-        (dti_filter, dti, [('out_bvec_file', 'bvec_file')]),
-        (dti_filter, dwidenoise, [('out_file', 'in_file')]),
+        (shells, dwimodel, [('out_data', 'bvals'),
+                            ('n_shells', 'n_shells')]),
+        (meta, dwimodel, [('out_bvec_file', 'bvec_file')]),
+        (drift, dwidenoise, [('out_full_file', 'in_file')]),
         (dmri_bmsk, dwidenoise, [('outputnode.out_mask', 'mask')]),
-        (dwidenoise, dti, [('out_file', 'in_file')]),
-        (dmri_bmsk, dti, [('outputnode.out_mask', 'brainmask')]),
+        (dwidenoise, dwimodel, [('out_file', 'in_file')]),
+        (dmri_bmsk, dwimodel, [('outputnode.out_mask', 'brain_mask')]),
         (meta, get_hmc_shells, [('out_bvec_file', 'in_bvec_file')]),
         (shells, get_hmc_shells, [('b_indices', 'indices')]),
         (hmcwf, get_hmc_shells, [('outputnode.out_file', 'in_file')]),
-        (dti, cc_mask, [('out_fa', 'in_fa'),
-                        ('out_cfa', 'in_cfa')]),
+        (dwimodel, cc_mask, [('out_fa', 'in_fa'),
+                             ('out_cfa', 'in_cfa')]),
         (meta, iqms_wf, [('qspace_neighbors', 'inputnode.qspace_neighbors')]),
         (averages, iqms_wf, [(('out_file', _first), 'inputnode.in_b0')]),
         (sp_mask, iqms_wf, [('out_mask', 'inputnode.spikes_mask')]),
         (piesno, iqms_wf, [('sigma', 'inputnode.piesno_sigma')]),
         (hmcwf, iqms_wf, [('outputnode.out_fd', 'inputnode.framewise_displacement')]),
-        (dti, iqms_wf, [('out_fa', 'inputnode.in_fa'),
-                        ('out_cfa', 'inputnode.in_cfa'),
-                        ('out_fa_nans', 'inputnode.in_fa_nans'),
-                        ('out_fa_degenerate', 'inputnode.in_fa_degenerate'),
-                        ('out_md', 'inputnode.in_md')]),
+        (dwimodel, iqms_wf, [('out_fa', 'inputnode.in_fa'),
+                             ('out_cfa', 'inputnode.in_cfa'),
+                             ('out_fa_nans', 'inputnode.in_fa_nans'),
+                             ('out_fa_degenerate', 'inputnode.in_fa_degenerate'),
+                             ('out_md', 'inputnode.in_md')]),
         (dmri_bmsk, iqms_wf, [('outputnode.out_mask', 'inputnode.brain_mask')]),
         (cc_mask, iqms_wf, [('out_mask', 'inputnode.cc_mask'),
                             ('wm_finalmask', 'inputnode.wm_mask')]),
@@ -281,8 +274,8 @@ def dmri_qc_workflow(name='dwiMRIQC'):
         (averages, dwi_report_wf, [('out_file', 'inputnode.in_avgmap')]),
         (stddev, dwi_report_wf, [('out_file', 'inputnode.in_stdmap')]),
         (drift, dwi_report_wf, [('out_full_file', 'inputnode.in_epi')]),
-        (dti, dwi_report_wf, [('out_fa', 'inputnode.in_fa'),
-                              ('out_md', 'inputnode.in_md')]),
+        (dwimodel, dwi_report_wf, [('out_fa', 'inputnode.in_fa'),
+                                   ('out_md', 'inputnode.in_md')]),
         (spatial_norm, dwi_report_wf, [('outputnode.epi_parc', 'inputnode.in_parcellation')]),
     ])
     # fmt: on

mriqc/workflows/diffusion/output.py

@@ -149,7 +149,7 @@ def init_dwi_report_wf(name='dwi_report_wf'):
     ds_report_md = pe.Node(
         DerivativesDataSink(
             base_directory=reportlets_dir,
-            desc='ad',
+            desc='md',
             datatype='figures',
         ),
         name='ds_report_md',