Merge pull request #1229 from nipreps/enh/dwidenoise-iqm

ENH: Noise floor estimated with PCA (``dwidenoise``) as an IQM

mriqc/interfaces/diffusion.py

@@ -104,6 +104,7 @@ class _DiffusionQCInputSpec(_BaseInterfaceInputSpec):
     wm_mask = File(exists=True, mandatory=True, desc='input probabilistic white-matter mask')
     cc_mask = File(exists=True, mandatory=True, desc='input binary mask of the corpus callosum')
     spikes_mask = File(exists=True, mandatory=True, desc='input binary mask of spiking voxels')
+    noise_floor = traits.Float(mandatory=True, desc='noise-floor map estimated by means of PCA')
     direction = traits.Enum(
         'all',
         'x',
@@ -142,6 +143,7 @@ class _DiffusionQCOutputSpec(TraitedSpec):
     fd = traits.Dict
     ndc = traits.Float
     sigma_cc = traits.Float
+    sigma_pca = traits.Float
     sigma_piesno = traits.Float
     spikes_ppm = traits.Dict
     # gsr = traits.Dict
@@ -223,13 +225,15 @@ def _run_interface(self, runtime):
         stats = aqc.summary_stats(b0data, rois)
         self._results['summary'] = stats
 
+        # CC mask SNR and std
         self._results['cc_snr'], cc_sigma = dqc.cc_snr(
             in_b0=b0data,
             dwi_shells=shelldata,
             cc_mask=ccdata,
             b_values=self.inputs.in_bval,
             b_vectors=self.inputs.in_bvec,
         )
+        self._results['sigma_cc'] = round(float(cc_sigma), 4)
 
         fa_nans_mask = np.asanyarray(nb.load(self.inputs.in_fa_nans).dataobj) > 0.0
         self._results['fa_nans'] = np.round(float(fa_nans_mask[mskdata > 0.5].mean()), 8) * 1e6
@@ -280,7 +284,9 @@ def _run_interface(self, runtime):
 
         # PIESNO
         self._results['sigma_piesno'] = round(self.inputs.piesno_sigma, 4)
-        self._results['sigma_cc'] = round(float(cc_sigma), 4)
+
+        # dwidenoise - Marchenko-Pastur PCA
+        self._results['sigma_pca'] = round(self.inputs.noise_floor, 4)
 
         self._results['out_qc'] = _flatten_dict(self._results)
 

mriqc/workflows/diffusion/base.py

@@ -271,9 +271,10 @@ def dmri_qc_workflow(name='dwiMRIQC'):
                            ('b_values', 'inputnode.b_values')]),
         (get_hmc_shells, iqms_wf, [('out_file', 'inputnode.in_shells'),
                                    ('out_bvec', 'inputnode.in_bvec')]),
+        (dwidenoise, iqms_wf, [('noise', 'inputnode.in_noise')]),
         (dwi_ref, spatial_norm, [('out_file', 'inputnode.epi_mean')]),
         (dmri_bmsk, spatial_norm, [('outputnode.out_mask', 'inputnode.epi_mask')]),
-        (dwidenoise, dwi_report_wf, [('noise', 'inputnode.in_noise')]),
+        (iqms_wf, dwi_report_wf, [('outputnode.noise_floor', 'inputnode.noise_floor')]),
         (shells, dwi_report_wf, [('b_dict', 'inputnode.in_bdict')]),
         (dmri_bmsk, dwi_report_wf, [('outputnode.out_mask', 'inputnode.brain_mask')]),
         (shells, dwi_report_wf, [('b_values', 'inputnode.in_shells')]),
@@ -324,6 +325,7 @@ def compute_iqms(name='ComputeIQMs'):
                 'in_fa_nans',
                 'in_fa_degenerate',
                 'in_md',
+                'in_noise',
                 'brain_mask',
                 'wm_mask',
                 'cc_mask',
@@ -340,11 +342,17 @@ def compute_iqms(name='ComputeIQMs'):
             fields=[
                 'out_file',
                 'meta_sidecar',
+                'noise_floor',
             ]
         ),
         name='outputnode',
     )
 
+    estimate_sigma = pe.Node(
+        niu.Function(function=_estimate_sigma),
+        name='estimate_sigma',
+    )
+
     meta = pe.Node(ReadSidecarJSON(index_db=config.execution.bids_database_dir), name='metadata')
 
     measures = pe.Node(DiffusionQC(), name='measures')
@@ -401,6 +409,10 @@ def compute_iqms(name='ComputeIQMs'):
         (datasink, outputnode, [('out_file', 'out_file')]),
         (meta, outputnode, [('out_dict', 'meta_sidecar')]),
         (measures, datasink, [('out_qc', 'root')]),
+        (inputnode, estimate_sigma, [('in_noise', 'in_file'),
+                                     ('brain_mask', 'mask')]),
+        (estimate_sigma, measures, [('out', 'noise_floor')]),
+        (estimate_sigma, outputnode, [('out', 'noise_floor')]),
     ])
     # fmt: on
     return workflow
@@ -639,3 +651,14 @@ def _all_but_first(inlist):
         return inlist[1:]
 
     return inlist
+
+
+def _estimate_sigma(in_file, mask):
+    import nibabel as nb
+    import numpy as np
+
+    msk = nb.load(mask).get_fdata() > 0.5
+    return round(
+        float(np.median(nb.load(in_file).get_fdata()[msk])),
+        6,
+    )

mriqc/workflows/diffusion/output.py

@@ -67,18 +67,13 @@ def init_dwi_report_wf(name='dwi_report_wf'):
                 'in_md',
                 'in_parcellation',
                 'in_bdict',
-                'in_noise',
+                'noise_floor',
                 'name_source',
             ]
         ),
         name='inputnode',
     )
 
-    estimate_sigma = pe.Node(
-        niu.Function(function=_estimate_sigma),
-        name='estimate_sigma',
-    )
-
     # Set FD threshold
     # inputnode.inputs.fd_thres = config.workflow.fd_thres
 
@@ -207,11 +202,9 @@ def _gen_entity(inlist):
         (inputnode, get_wm, [('in_parcellation', 'in_file')]),
         (inputnode, plot_heatmap, [('in_epi', 'in_file'),
                                    ('in_fa', 'scalarmap'),
-                                   ('in_bdict', 'b_indices')]),
+                                   ('in_bdict', 'b_indices'),
+                                   ('noise_floor', 'sigma')]),
         (inputnode, ds_report_hm, [('name_source', 'source_file')]),
-        (inputnode, estimate_sigma, [('in_noise', 'in_file'),
-                                     ('brain_mask', 'mask')]),
-        (estimate_sigma, plot_heatmap, [('out', 'sigma')]),
         (get_wm, plot_heatmap, [('out', 'mask_file')]),
         (plot_heatmap, ds_report_hm, [('out_file', 'in_file')]),
 
@@ -425,14 +418,3 @@ def _get_wm(in_file, radius=2):
         hdr,
     ).to_filename(out_wm)
     return out_wm
-
-
-def _estimate_sigma(in_file, mask):
-    import nibabel as nb
-    import numpy as np
-
-    msk = np.asanyarray(nb.load(mask).dataobj) > 0.5
-
-    return float(
-        np.median(nb.load(in_file).get_fdata()[msk])
-    )