Merge pull request #1115 from nipreps/enh/dwi-continuation

ENH: Add a basic DTI fitting into the diffusion workflow

mriqc/data/bootstrap-dwi.yml

@@ -31,6 +31,12 @@ sections:
 - name: Summary
   reportlets:
   - bids: {datatype: figures, desc: summary, extension: [.html]}
+  - bids: {datatype: figures, desc: fa}
+    caption: Reconstructed FA map using Dipy's free-water DTI model.
+    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
 - name: DWI shells
   ordering: bval
   reportlets:

mriqc/interfaces/diffusion.py

@@ -346,6 +346,96 @@ def _run_interface(self, runtime):
         return runtime
 
 
+class _DipyDTIInputSpec(_BaseInterfaceInputSpec):
+    in_file = File(exists=True, mandatory=True, desc="dwi file")
+    brainmask = File(exists=True, desc="brain mask file")
+    bvals = traits.List(traits.Float, mandatory=True, desc="bval table")
+    bvec_file = File(exists=True, mandatory=True, desc="b-vectors")
+    free_water_model = traits.Bool(False, usedefault=True, desc="use free water model")
+
+
+class _DipyDTIOutputSpec(_TraitedSpec):
+    out_fa = File(exists=True, desc="output FA file")
+    out_md = File(exists=True, desc="output MD file")
+
+
+class DipyDTI(SimpleInterface):
+    """Split a DWI dataset into ."""
+
+    input_spec = _DipyDTIInputSpec
+    output_spec = _DipyDTIOutputSpec
+
+    def _run_interface(self, runtime):
+        from dipy.core.gradients import gradient_table_from_bvals_bvecs
+        from dipy.reconst.dti import TensorModel
+        from dipy.reconst.fwdti import FreeWaterTensorModel
+        from nipype.utils.filemanip import fname_presuffix
+
+        gtab = gradient_table_from_bvals_bvecs(
+            bvals=self.inputs.bvals,
+            bvecs=np.loadtxt(self.inputs.bvec_file).T,
+        )
+
+        img = nb.load(self.inputs.in_file)
+        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
+
+        DTIModel = FreeWaterTensorModel if self.inputs.free_water_model else TensorModel
+
+        # Fit DIT
+        fwdtifit = DTIModel(gtab).fit(
+            data,
+            mask=brainmask,
+        )
+
+        # Extract the FA
+        fa_data = np.array(fwdtifit.fa, dtype="float32")
+        fa_data[np.isnan(fa_data)] = 0
+        fa_data = np.clip(fa_data, 0, 1)
+
+        fa_nii = nb.Nifti1Image(
+            fa_data,
+            img.affine,
+            None,
+        )
+
+        fa_nii.header.set_xyzt_units("mm")
+        fa_nii.header.set_intent("estimate", name="Fractional Anisotropy (FA)")
+        fa_nii.header["cal_max"] = 1.0
+        fa_nii.header["cal_min"] = 0.0
+
+        self._results["out_fa"] = fname_presuffix(
+            self.inputs.in_file,
+            suffix="fa",
+            newpath=runtime.cwd,
+        )
+
+        fa_nii.to_filename(self._results["out_fa"])
+
+        # Extract the AD
+        self._results["out_md"] = fname_presuffix(
+            self.inputs.in_file,
+            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)")
+        nb.Nifti1Image(
+            ad_data,
+            img.affine,
+            ad_hdr
+        ).to_filename(self._results["out_md"])
+
+        return runtime
+
+
 def _rms(estimator, X):
     """
     Callable to pass to GridSearchCV that will calculate a distance score.

mriqc/workflows/diffusion/base.py

@@ -67,6 +67,7 @@ def dmri_qc_workflow(name="dwiMRIQC"):
     from niworkflows.workflows.epi.refmap import init_epi_reference_wf
     from mriqc.interfaces.diffusion import (
         CorrectSignalDrift,
+        DipyDTI,
         ExtractB0,
         NumberOfShells,
         ReadDWIMetadata,
@@ -161,13 +162,19 @@ def dmri_qc_workflow(name="dwiMRIQC"):
         iterfield=["in_weights"],
     )
 
-    # 6. EPI to MNI registration
+    # 6. Fit DTI model
+    dti = pe.Node(
+        DipyDTI(free_water_model=False),
+        name="dti",
+    )
+
+    # 7. EPI to MNI registration
     ema = epi_mni_align()
 
-    # 7. Compute IQMs
+    # 8. Compute IQMs
     iqmswf = compute_iqms()
 
-    # 8. Generate outputs
+    # 9. Generate outputs
     dwi_report_wf = init_dwi_report_wf()
 
     # fmt: off
@@ -201,13 +208,19 @@ def dmri_qc_workflow(name="dwiMRIQC"):
         (drift, stddev, [("out_full_file", "in_file")]),
         (shells, averages, [("b_masks", "in_weights")]),
         (shells, stddev, [("b_masks", "in_weights")]),
+        (shells, dti, [("out_data", "bvals")]),
+        (meta, dti, [("out_bvec_file", "bvec_file")]),
+        (drift, dti, [("out_full_file", "in_file")]),
+        (dmri_bmsk, dti, [("outputnode.out_mask", "brainmask")]),
         (hmcwf, outputnode, [("outputnode.out_fd", "out_fd")]),
         (shells, iqmswf, [("n_shells", "inputnode.n_shells"),
                           ("b_values", "inputnode.b_values")]),
         (dmri_bmsk, dwi_report_wf, [("outputnode.out_mask", "inputnode.brainmask")]),
         (shells, dwi_report_wf, [("b_values", "inputnode.in_shells")]),
         (averages, dwi_report_wf, [("out_file", "inputnode.in_avgmap")]),
         (stddev, dwi_report_wf, [("out_file", "inputnode.in_stdmap")]),
+        (dti, dwi_report_wf, [("out_fa", "inputnode.in_fa"),
+                              ("out_md", "inputnode.in_md")]),
     ])
     # fmt: on
     return workflow

mriqc/workflows/diffusion/output.py

@@ -60,6 +60,8 @@ def init_dwi_report_wf(name="dwi_report_wf"):
                 "in_avgmap",
                 "in_stdmap",
                 "in_shells",
+                "in_fa",
+                "in_md",
                 # "in_ras",
                 # "hmc_epi",
                 # "epi_mean",
@@ -84,6 +86,15 @@ def init_dwi_report_wf(name="dwi_report_wf"):
     # Set FD threshold
     # inputnode.inputs.fd_thres = config.workflow.fd_thres
 
+    mosaic_fa = pe.Node(
+        PlotMosaic(cmap="Greys_r"),
+        name="mosaic_fa",
+    )
+    mosaic_md = pe.Node(
+        PlotMosaic(cmap="Greys_r"),
+        name="mosaic_md",
+    )
+
     mosaic_mean = pe.MapNode(
         PlotMosaic(cmap="Greys_r"),
         name="mosaic_mean",
@@ -108,8 +119,9 @@ def init_dwi_report_wf(name="dwi_report_wf"):
     if config.workflow.species.lower() in ("rat", "mouse"):
         mosaic_mean.inputs.view = ["coronal", "axial"]
         mosaic_stddev.inputs.view = ["coronal", "axial"]
-        # mosaic_zoom.inputs.view = ["coronal", "axial"]
         mosaic_noise.inputs.view = ["coronal", "axial"]
+        mosaic_fa.inputs.view = ["coronal", "axial"]
+        mosaic_md.inputs.view = ["coronal", "axial"]
 
     ds_report_mean = pe.MapNode(
         DerivativesDataSink(
@@ -147,6 +159,26 @@ def init_dwi_report_wf(name="dwi_report_wf"):
         iterfield=["in_file", "bval"],
     )
 
+    ds_report_fa = pe.Node(
+        DerivativesDataSink(
+            base_directory=reportlets_dir,
+            desc="fa",
+            datatype="figures",
+        ),
+        name="ds_report_fa",
+        run_without_submitting=True,
+    )
+
+    ds_report_md = pe.Node(
+        DerivativesDataSink(
+            base_directory=reportlets_dir,
+            desc="ad",
+            datatype="figures",
+        ),
+        name="ds_report_md",
+        run_without_submitting=True,
+    )
+
     def _gen_entity(inlist):
         return ["00000"] + [f"{int(round(bval, 0)):05d}" for bval in inlist]
 
@@ -157,15 +189,23 @@ def _gen_entity(inlist):
         (inputnode, mosaic_stddev, [("in_stdmap", "in_file"),
                                     ("brainmask", "bbox_mask_file")]),
         (inputnode, mosaic_noise, [("in_avgmap", "in_file")]),
+        (inputnode, mosaic_fa, [("in_fa", "in_file"),
+                                ("brainmask", "bbox_mask_file")]),
+        (inputnode, mosaic_md, [("in_md", "in_file"),
+                                ("brainmask", "bbox_mask_file")]),
         (inputnode, ds_report_mean, [("name_source", "source_file"),
                                      (("in_shells", _gen_entity), "bval")]),
         (inputnode, ds_report_stdev, [("name_source", "source_file"),
                                       (("in_shells", _gen_entity), "bval")]),
         (inputnode, ds_report_noise, [("name_source", "source_file"),
                                       (("in_shells", _gen_entity), "bval")]),
+        (inputnode, ds_report_fa, [("name_source", "source_file")]),
+        (inputnode, ds_report_md, [("name_source", "source_file")]),
         (mosaic_mean, ds_report_mean, [("out_file", "in_file")]),
         (mosaic_stddev, ds_report_stdev, [("out_file", "in_file")]),
         (mosaic_noise, ds_report_noise, [("out_file", "in_file")]),
+        (mosaic_fa, ds_report_fa, [("out_file", "in_file")]),
+        (mosaic_md, ds_report_md, [("out_file", "in_file")]),
     ])
     # fmt: on