From f1af381cb0562ee5ff9f15295bbab32a400efac2 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Wed, 3 Apr 2024 11:30:32 +0200 Subject: [PATCH] ENH: Noise floor estimated with PCA (``dwidenoise``) as an IQM --- mriqc/interfaces/diffusion.py | 8 +++++++- mriqc/workflows/diffusion/base.py | 25 ++++++++++++++++++++++++- mriqc/workflows/diffusion/output.py | 24 +++--------------------- 3 files changed, 34 insertions(+), 23 deletions(-) diff --git a/mriqc/interfaces/diffusion.py b/mriqc/interfaces/diffusion.py index f204058e..3ab5de3b 100644 --- a/mriqc/interfaces/diffusion.py +++ b/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,6 +225,7 @@ 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, @@ -230,6 +233,7 @@ def _run_interface(self, runtime): 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) diff --git a/mriqc/workflows/diffusion/base.py b/mriqc/workflows/diffusion/base.py index 609a2564..6dc008d5 100644 --- a/mriqc/workflows/diffusion/base.py +++ b/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, + ) diff --git a/mriqc/workflows/diffusion/output.py b/mriqc/workflows/diffusion/output.py index de4062c4..e4b9e2af 100644 --- a/mriqc/workflows/diffusion/output.py +++ b/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]) - )