diff --git a/mriqc/data/bootstrap-dwi.yml b/mriqc/data/bootstrap-dwi.yml index a764f68e..999cd11b 100644 --- a/mriqc/data/bootstrap-dwi.yml +++ b/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: diff --git a/mriqc/interfaces/diffusion.py b/mriqc/interfaces/diffusion.py index 3ab5de3b..e44da956 100644 --- a/mriqc/interfaces/diffusion.py +++ b/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 diff --git a/mriqc/workflows/diffusion/base.py b/mriqc/workflows/diffusion/base.py index 6dc008d5..2cfef125 100644 --- a/mriqc/workflows/diffusion/base.py +++ b/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 diff --git a/mriqc/workflows/diffusion/output.py b/mriqc/workflows/diffusion/output.py index e4b9e2af..0fe963a7 100644 --- a/mriqc/workflows/diffusion/output.py +++ b/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',