Merge pull request #1119 from nipreps/enh/add-heatmaps

ENH: Incorporate new NiReports' DWI heatmaps

mriqc/data/bootstrap-dwi.yml

@@ -31,6 +31,12 @@ sections:
 - name: Summary
   reportlets:
   - bids: {datatype: figures, desc: summary, extension: [.html]}
+  - bids: {datatype: figures, desc: heatmap}
+    caption: This visualization divides the data by shells, and shows the joint distribution
+      of SNR vs. FA. At the bottom, the distributions are marginalized for SNR.
+      Please note that the figures of SNR provided are calculated with a coarse estimation of
+      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.
     subtitle: Fractional anisotropy (FA) map

mriqc/interfaces/diffusion.py

@@ -134,6 +134,9 @@ class _NumberOfShellsOutputSpec(_TraitedSpec):
         traits.List(traits.Int, minlen=1),
         minlen=1,
         desc="b-value-wise masks")
+    b_dict = traits.Dict(
+        traits.Int, traits.List(traits.Int), desc="b-values dictionary"
+    )
 
 
 class NumberOfShells(SimpleInterface):
@@ -200,6 +203,11 @@ def _run_interface(self, runtime):
             np.atleast_1d(np.squeeze(np.argwhere(b_mask)).astype(int)).tolist()
             for b_mask in self._results["b_masks"]
         ]
+
+        self._results["b_dict"] = {
+            int(round(k, 0)): value
+            for k, value in zip([0] + self._results["b_values"], self._results["b_indices"])
+        }
         return runtime
 
 
@@ -346,12 +354,74 @@ def _run_interface(self, runtime):
         return runtime
 
 
+class _FilterShellsInputSpec(_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")
+    b_threshold = traits.Float(1100, usedefault=True, desc="b-values threshold")
+
+
+class _FilterShellsOutputSpec(_TraitedSpec):
+    out_file = File(exists=True, desc="filtered DWI file")
+    out_bvals = traits.List(traits.Float, desc="filtered bvalues")
+    out_bvec_file = File(exists=True, desc="filtered bvecs file")
+    out_bval_file = File(exists=True, desc="filtered bvals file")
+
+
+class FilterShells(SimpleInterface):
+    """Extract DWIs below a given b-value threshold."""
+
+    input_spec = _FilterShellsInputSpec
+    output_spec = _FilterShellsOutputSpec
+
+    def _run_interface(self, runtime):
+        from nipype.utils.filemanip import fname_presuffix
+
+        bvals = np.array(self.inputs.bvals)
+        bval_mask = bvals < self.inputs.b_threshold
+        bvecs = np.loadtxt(self.inputs.bvec_file)[:, bval_mask]
+
+        self._results["out_bvals"] = bvals[bval_mask].astype(float).tolist()
+        self._results["out_bvec_file"] = fname_presuffix(
+            self.inputs.in_file,
+            suffix="_dti.bvec",
+            newpath=runtime.cwd,
+            use_ext=False,
+        )
+        np.savetxt(self._results["out_bvec_file"], bvecs)
+
+        self._results["out_bval_file"] = fname_presuffix(
+            self.inputs.in_file,
+            suffix="_dti.bval",
+            newpath=runtime.cwd,
+            use_ext=False,
+        )
+        np.savetxt(self._results["out_bval_file"], bvals)
+
+        self._results["out_file"] = fname_presuffix(
+            self.inputs.in_file,
+            suffix="_dti",
+            newpath=runtime.cwd,
+        )
+
+        dwi_img = nb.load(self.inputs.in_file)
+        data = np.array(dwi_img.dataobj, dtype=dwi_img.header.get_data_dtype())[..., bval_mask]
+        dwi_img.__class__(
+            data,
+            dwi_img.affine,
+            dwi_img.header,
+        ).to_filename(self._results["out_file"])
+
+        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")
+    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")
 
 
 class _DipyDTIOutputSpec(_TraitedSpec):
@@ -371,13 +441,16 @@ def _run_interface(self, runtime):
         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=self.inputs.bvals,
-            bvecs=np.loadtxt(self.inputs.bvec_file).T,
+            bvals=bvals[bval_mask],
+            bvecs=np.loadtxt(self.inputs.bvec_file).T[bval_mask],
         )
 
         img = nb.load(self.inputs.in_file)
-        data = img.get_fdata(dtype="float32")
+        data = img.get_fdata(dtype="float32")[..., bval_mask]
 
         brainmask = np.ones_like(data[..., 0], dtype=bool)
 

mriqc/workflows/diffusion/base.py

@@ -63,12 +63,14 @@ def dmri_qc_workflow(name="dwiMRIQC"):
 
     """
     from nipype.interfaces.afni import Volreg
+    from nipype.interfaces.mrtrix3.preprocess import DWIDenoise
     from niworkflows.interfaces.header import SanitizeImage
     from niworkflows.workflows.epi.refmap import init_epi_reference_wf
     from mriqc.interfaces.diffusion import (
         CorrectSignalDrift,
         DipyDTI,
         ExtractB0,
+        FilterShells,
         NumberOfShells,
         ReadDWIMetadata,
         WeightedStat,
@@ -163,6 +165,15 @@ def dmri_qc_workflow(name="dwiMRIQC"):
     )
 
     # 6. Fit DTI model
+    dti_filter = pe.Node(FilterShells(), name="dti_filter")
+    dwidenoise = pe.Node(
+        DWIDenoise(
+            noise="noisemap.nii.gz",
+            nthreads=config.nipype.omp_nthreads,
+        ),
+        name="dwidenoise",
+        nprocs=config.nipype.omp_nthreads,
+    )
     dti = pe.Node(
         DipyDTI(free_water_model=False),
         name="dti",
@@ -208,19 +219,28 @@ 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")]),
+        (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")]),
+        (dmri_bmsk, dwidenoise, [("outputnode.out_mask", "mask")]),
+        (dwidenoise, dti, [("out_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")]),
+        (dwidenoise, dwi_report_wf, [("noise", "inputnode.in_noise")]),
+        (shells, dwi_report_wf, [("b_dict", "inputnode.in_bdict")]),
         (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")]),
+        (drift, dwi_report_wf, [("out_full_file", "inputnode.in_epi")]),
         (dti, dwi_report_wf, [("out_fa", "inputnode.in_fa"),
                               ("out_md", "inputnode.in_md")]),
+        (ema, dwi_report_wf, [("outputnode.epi_parc", "inputnode.in_parcellation")]),
     ])
     # fmt: on
     return workflow