diff --git a/mriqc/interfaces/diffusion.py b/mriqc/interfaces/diffusion.py index 91aae38b..4ee97d46 100644 --- a/mriqc/interfaces/diffusion.py +++ b/mriqc/interfaces/diffusion.py @@ -21,18 +21,24 @@ # https://www.nipreps.org/community/licensing/ # """Interfaces for manipulating DWI data.""" +from __future__ import annotations + import nibabel as nb import numpy as np +from dipy.reconst.dti import TensorFit from nipype.interfaces.base import ( - BaseInterfaceInputSpec as _BaseInterfaceInputSpec, -) -from nipype.interfaces.base import ( + BaseInterfaceInputSpec, File, + InputMultiObject, OutputMultiObject, SimpleInterface, + TraitedSpec, isdefined, traits, ) +from nipype.interfaces.base import ( + BaseInterfaceInputSpec as _BaseInterfaceInputSpec, +) from nipype.interfaces.base import ( TraitedSpec as _TraitedSpec, ) @@ -428,6 +434,8 @@ class _DipyDTIInputSpec(_BaseInterfaceInputSpec): class _DipyDTIOutputSpec(_TraitedSpec): out_fa = File(exists=True, desc='output FA file') + out_fa_clean = File(exists=True, desc='output FA file after fixing degenerate vectors') + out_cfa = File(exists=True, desc='output color FA file') out_md = File(exists=True, desc='output MD file') @@ -440,6 +448,7 @@ class DipyDTI(SimpleInterface): def _run_interface(self, runtime): from dipy.core.gradients import gradient_table_from_bvals_bvecs from dipy.reconst.dti import TensorModel + from dipy.reconst.dti import fractional_anisotropy, color_fa from dipy.reconst.fwdti import FreeWaterTensorModel from nipype.utils.filemanip import fname_presuffix @@ -468,9 +477,8 @@ def _run_interface(self, runtime): ) # 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_data = fractional_anisotropy(fwdtifit.evals) + fa_data[np.isnan(fa_data)] = -1000 # Arbitrary value to mark up NaNs fa_nii = nb.Nifti1Image( fa_data, @@ -480,8 +488,6 @@ def _run_interface(self, runtime): 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, @@ -491,6 +497,49 @@ def _run_interface(self, runtime): fa_nii.to_filename(self._results['out_fa']) + # Clamp the FA to remove degenerate tensors and round for stability + fa_clean = np.clip(np.round(fa_data, 3), 0, 1) + + fa_clean_nii = nb.Nifti1Image( + fa_clean, + img.affine, + None, + ) + + fa_clean_nii.header.set_xyzt_units('mm') + fa_clean_nii.header.set_intent( + 'estimate', + name='Fractional Anisotropy (FA) after removing degenerate vectors') + fa_clean_nii.header['cal_max'] = 1.0 + fa_clean_nii.header['cal_min'] = 0.0 + + self._results['out_fa_clean'] = fname_presuffix( + self.inputs.in_file, + suffix='desc-clean_fa', + newpath=runtime.cwd, + ) + fa_clean_nii.to_filename(self._results['out_fa_clean']) + + # Extract the color FA + cfa_data = color_fa(fa_clean, fwdtifit.evecs) + cfa_nii = nb.Nifti1Image( + cfa_data, + img.affine, + None, + ) + + cfa_nii.header.set_xyzt_units('mm') + cfa_nii.header.set_intent('estimate', name='Fractional Anisotropy (FA)') + cfa_nii.header['cal_max'] = 1.0 + cfa_nii.header['cal_min'] = 0.0 + + self._results['out_cfa'] = fname_presuffix( + self.inputs.in_file, + suffix='cfa', + newpath=runtime.cwd, + ) + cfa_nii.to_filename(self._results['out_cfa']) + # Extract the AD self._results['out_md'] = fname_presuffix( self.inputs.in_file, @@ -511,6 +560,110 @@ def _run_interface(self, runtime): return runtime +class _DiffusionQCInputSpec(BaseInterfaceInputSpec): + in_b0 = File(exists=True, mandatory=True, desc='input b=0 average') + in_shells = InputMultiObject( + File(exists=True), + mandatory=True, + desc='DWI data after HMC and split by shells (indexed by in_bval)' + ) + in_bvec = File(exists=True, mandatory=True, desc='input motion corrected file') + in_bval = traits.List(traits.Float, minlen=1, desc='b-values') + in_bvec = traits.List( + traits.Tuple(traits.Float, traits.Float, traits.Float), + minlen=1, + desc='b-vectors', + ) + in_mask = File(exists=True, mandatory=True, desc='input probabilistic brain mask') + in_fa = File(exists=True, mandatory=True, desc='input FA map') + in_md = File(exists=True, mandatory=True, desc='input MD map') + direction = traits.Enum( + 'all', + 'x', + 'y', + '-x', + '-y', + usedefault=True, + desc='direction for GSR computation', + ) + in_fwhm = traits.List(traits.Float, mandatory=True, desc='smoothness estimated with AFNI') + + +class _DiffusionQCOutputSpec(TraitedSpec): + fber = traits.Float + efc = traits.Float + snr = traits.Float + gsr = traits.Dict + tsnr = traits.Float + fd = traits.Dict + fwhm = traits.Dict(desc='full width half-maximum measure') + size = traits.Dict + spacing = traits.Dict + summary = traits.Dict + + out_qc = traits.Dict(desc='output flattened dictionary with all measures') + + +class DiffusionQC(SimpleInterface): + """Computes :abbr:`QC (Quality Control)` measures on the input DWI EPI scan.""" + + input_spec = _DiffusionQCInputSpec + output_spec = _DiffusionQCOutputSpec + + def _run_interface(self, runtime): + from mriqc.qc import anatomical as aqc + from mriqc.qc import diffusion as dqc + + # Get the mean EPI data and get it ready + b0nii = nb.load(self.inputs.in_b0) + b0data = np.round( + np.nan_to_num(np.asanyarray(b0nii.dataobj)), + 2, + ) + b0data[b0data < 0] = 0 + + # Get the FA data and get it ready + fanii = nb.load(self.inputs.in_fa) + fadata = np.round( + np.nan_to_num(np.asanyarray(fanii.dataobj)), + 2, + ) + + # Get EPI data (with hmc done) and get it ready + hmcnii = nb.load(self.inputs.in_hmc) + hmcdata = np.round( + np.nan_to_num(np.asanyarray(hmcnii.dataobj)), + 2, + ) + hmcdata[hmcdata < 0] = 0 + + # Get brain mask data + msknii = nb.load(self.inputs.in_mask) + mskdata = np.round( # Protect the thresholding with a rounding for stability + np.asanyarray(msknii.dataobj), + 1, + ) > 0 + if np.sum(mskdata) < 100: + raise RuntimeError( + 'Detected less than 100 voxels belonging to the brain mask. ' + 'MRIQC failed to process this dataset.' + ) + + # Summary stats + rois = { + 'fg': mskdata, + 'bg': 1.0 - mskdata, + } + stats = aqc.summary_stats(b0data, rois) + self._results['summary'] = stats + + self._results['cc_snr'] = dqc.cc_snr( + fadata + ) + + return runtime + + def _rms(estimator, X): """ Callable to pass to GridSearchCV that will calculate a distance score. @@ -545,3 +698,154 @@ def _extract_b0(in_file, b0_ixs, out_path=None): def _exp_func(t, A, K, C): return A * np.exp(K * t) + C + + +def get_spike_mask( + data: np.ndarray, + z_threshold: float = 3.0, + grouping_vals: np.ndarray | None = None, + bmag: int | None = None, +) -> np.ndarray: + """ + Creates a binary mask classifying voxels in the data array as spike or non-spike. + + This function identifies voxels with signal intensities exceeding a threshold based + on standard deviations above the mean. The threshold can be applied globally to + the entire data array, or it can be calculated for groups of voxels defined by + the `grouping_vals` parameter. + + Parameters + ---------- + data : ndarray (float, 3D or 4D) + The data array to be thresholded. + z_threshold : float, optional (default=3.0) + The number of standard deviations to use above the mean as the threshold + multiplier. + grouping_vals : ndarray (int, same shape as data or 3D), optional + If provided, this array is used to group voxels for thresholding. Voxels + with the same value in `grouping_vals` are considered to belong to the same + group. The threshold will be calculated independently for each group. + - If `grouping_vals` has the same shape as `data` (4D), it is assumed to be + a mask where each voxel value indicates the group it belongs to. + - If `grouping_vals` has a 3D shape, it is assumed to represent b-values + corresponding to each voxel in the 4D `data` array. In this case, voxels + with the same b-value are grouped together. + bmag : int, optional + The order of magnitude for b-value rounding (used only if + `grouping_vals` is provided as b-values). Default: None (derived from max b-value). + + Returns: + ------- + ndarray (bool, same shape as data) + A binary mask where True values indicate voxels classified as spikes and + False values indicate non-spikes. The mask has the same shape as the input + data array. + + """ + from dipy.core.gradients import round_bvals, unique_bvals_magnitude + + if grouping_vals is None: + threshold = np.round((z_threshold * np.std(data)) + np.mean(data), 3) + spike_mask = np.round(data, 3) > threshold + return spike_mask + + threshold_mask = np.zeros(data.shape) + + rounded_grouping_vals = round_bvals(grouping_vals, bmag) + gvals = unique_bvals_magnitude(grouping_vals, bmag) + + if grouping_vals.shape == data.shape: + for gval in gvals: + gval_data = data[rounded_grouping_vals == gval] + gval_threshold = ((z_threshold * np.std(gval_data)) + + np.mean(gval_data)) + threshold_mask[rounded_grouping_vals == gval] = ( + gval_threshold * np.ones(gval_data.shape)) + else: + for gval in gvals: + gval_data = data[..., rounded_grouping_vals == gval] + gval_threshold = ((z_threshold * np.std(gval_data)) + + np.mean(gval_data)) + threshold_mask[..., rounded_grouping_vals == gval] = ( + gval_threshold * np.ones(gval_data.shape)) + + spike_mask = data > threshold_mask + + return spike_mask + + +def noise_b0(data, gtab, mask=None): + """ + Estimate noise in raw dMRI based on b0 variance. + + Parameters + ---------- + """ + if mask is None: + mask = np.ones(data.shape[:3], dtype=bool) + b0 = data[..., ~gtab.b0s_mask] + return np.percentile(np.var(b0[mask], -1), (25, 50, 75)) + + +def segment_corpus_callosum( + in_cfa: np.ndarray, + mask: np.ndarray, + threshold: tuple[float, float, float, float, float, float] = ( + 0.6, + 1, + 0, + 0.1, + 0, + 0.1, + ), +) -> tuple[np.ndarray, np.ndarray]: + """ + Segments the corpus callosum (CC) from a color FA map. + + Parameters + ---------- + in_cfa : :obj:`~numpy.ndarray` + The color FA (cFA) map. + mask : :obj:`~numpy.ndarray` (bool, 3D) + A brain mask used to define the initial bounding box. + threshold : :obj:`tuple`, optional + An iterable that defines the minimum and maximum values to use for + the thresholding of the cFA. Values are specified as + (R_min, R_max, G_min, G_max, B_min, B_max). + + Returns + ------- + cc_mask: :obj:`~numpy.ndarray` + The final binary mask of the segmented CC. + + Notes + ----- + This implementation was derived from + :obj:`dipy.segment.mask.segment_from_cfa`. + The CC mask is then cleaned-up for spurious off voxels with + :obj:`dipy.segment.mask.clean_cc_mask` + + """ + from dipy.segment.mask import bounding_box, clean_cc_mask + + # Prepare a bounding box of the CC + cc_box = np.zeros_like(mask, dtype=bool) + mins, maxs = bounding_box(mask) # mask needs to be volume + mins = np.array(mins) + maxs = np.array(maxs) + diff = (maxs - mins) // 4 + bounds_min = mins + diff + bounds_max = maxs - diff + cc_box[ + bounds_min[0]:bounds_max[0], + bounds_min[1]:bounds_max[1], + bounds_min[2]:bounds_max[2] + ] = True + + include = ( + (in_cfa >= threshold[0::2]) + & (in_cfa <= threshold[1::2]) + & cc_box[..., None] + ) + cc_mask = clean_cc_mask(np.all(include, axis=-1)) + return cc_mask