From e346e1ccba846d09ab29a23d5366cc9882fc572f Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Wed, 3 Apr 2024 17:12:09 +0200 Subject: [PATCH] enh: add b-vector angular deviations as IQMs References: #1216. --- mriqc/interfaces/diffusion.py | 96 ++++++++++++++++++++++++++++++- mriqc/workflows/diffusion/base.py | 44 ++++++++------ 2 files changed, 122 insertions(+), 18 deletions(-) diff --git a/mriqc/interfaces/diffusion.py b/mriqc/interfaces/diffusion.py index e44da956..b0cb0567 100644 --- a/mriqc/interfaces/diffusion.py +++ b/mriqc/interfaces/diffusion.py @@ -51,13 +51,14 @@ __all__ = ( 'CCSegmentation', 'CorrectSignalDrift', - 'DiffusionQC', 'DiffusionModel', + 'DiffusionQC', 'ExtractOrientations', 'FilterShells', 'NumberOfShells', 'PIESNO', 'ReadDWIMetadata', + 'RotateVectors', 'SpikingVoxelsMask', 'SplitShells', 'WeightedStat', @@ -90,6 +91,18 @@ class _DiffusionQCInputSpec(_BaseInterfaceInputSpec): minlen=1, desc='a list of shell-wise splits of b-vectors lists -- first list are b=0', ) + in_bvec_rotated = traits.List( + traits.Tuple(traits.Float, traits.Float, traits.Float), + mandatory=True, + minlen=1, + desc='b-vectors after rotating by the head-motion correction transform', + ) + in_bvec_diff = traits.List( + traits.Float, + mandatory=True, + minlen=1, + desc='list of angle deviations from the original b-vectors table', + ) in_fa = File(exists=True, mandatory=True, desc='input FA map') in_fa_nans = File(exists=True, mandatory=True, desc='binary mask of NaN values in the "raw" FA map') @@ -135,6 +148,7 @@ class _DiffusionQCInputSpec(_BaseInterfaceInputSpec): class _DiffusionQCOutputSpec(TraitedSpec): + bdiffs = traits.Dict cc_snr = traits.Dict efc = traits.Dict fa_degenerate = traits.Float @@ -288,8 +302,16 @@ def _run_interface(self, runtime): # dwidenoise - Marchenko-Pastur PCA self._results['sigma_pca'] = round(self.inputs.noise_floor, 4) - self._results['out_qc'] = _flatten_dict(self._results) + # rotated b-vecs deviations + diffs = np.array(self.inputs.in_bvec_diff) + self._results['bdiffs'] = { + 'mean': round(float(diffs[diffs > 1e-4].mean()), 4), + 'median': round(float(np.median(diffs[diffs > 1e-4])), 4), + 'max': round(float(diffs[diffs > 1e-4].max()), 4), + 'min': round(float(diffs[diffs > 1e-4].min()), 4), + } + self._results['out_qc'] = _flatten_dict(self._results) return runtime @@ -1093,6 +1115,76 @@ def _run_interface(self, runtime): return runtime +class _RotateVectorsInputSpec(_BaseInterfaceInputSpec): + in_file = File( + exists=True, + mandatory=True, + desc='TSV file containing original b-vectors and b-values', + ) + reference = File( + exists=True, + mandatory=True, + desc='dwi-related file providing the reference affine', + ) + transforms = File(exists=True, desc='list of head-motion transforms') + + +class _RotateVectorsOutputSpec(_TraitedSpec): + out_bvec = traits.List( + traits.Tuple(traits.Float, traits.Float, traits.Float), + minlen=1, + desc='rotated b-vectors', + ) + out_diff = traits.List( + traits.Float, + minlen=1, + desc='angles in radians between new b-vectors and the original ones', + ) + + +class RotateVectors(SimpleInterface): + """Extract all b=0 volumes from a dwi series.""" + + input_spec = _RotateVectorsInputSpec + output_spec = _RotateVectorsOutputSpec + + def _run_interface(self, runtime): + from nitransforms.linear import load + + vox2ras = nb.load(self.inputs.reference).affine + ras2vox = np.linalg.inv(vox2ras) + + ijk = np.loadtxt(self.inputs.in_file).T + nonzero = np.linalg.norm(ijk, axis=1) > 1e-3 + + xyz = (vox2ras[:3, :3] @ ijk.T).T + + # Unit vectors in RAS coordinates + xyz_norms = np.linalg.norm(xyz, axis=1) + xyz[nonzero] = xyz[nonzero] / xyz_norms[nonzero, np.newaxis] + + hmc_rot = load(self.inputs.transforms).matrix[:, :3, :3] + ijk_rotated = ( + ras2vox[:3, :3] + @ np.einsum('ijk,ik->ij', hmc_rot, xyz).T + ).T.astype('float32') + ijk_rotated_norm = np.linalg.norm(ijk_rotated, axis=1) + ijk_rotated[nonzero] = ijk_rotated[nonzero] / ijk_rotated_norm[nonzero, np.newaxis] + ijk_rotated[~nonzero] = ijk[~nonzero] + + self._results['out_bvec'] = list( + zip(ijk_rotated[:, 0], ijk_rotated[:, 1], ijk_rotated[:, 2]) + ) + + diffs = np.zeros_like(ijk[:, 0]) + diffs[nonzero] = np.arccos( + np.clip(np.einsum('ij, ij->i', ijk[nonzero], ijk_rotated[nonzero]), -1.0, 1.0) + ) + self._results['out_diff'] = [round(float(v), 6) for v in diffs] + + return runtime + + def _rms(estimator, X): """ Callable to pass to GridSearchCV that will calculate a distance score. diff --git a/mriqc/workflows/diffusion/base.py b/mriqc/workflows/diffusion/base.py index 2cfef125..f5afb48b 100644 --- a/mriqc/workflows/diffusion/base.py +++ b/mriqc/workflows/diffusion/base.py @@ -110,11 +110,6 @@ def dmri_qc_workflow(name='dwiMRIQC'): name='datalad_get', ) - # outputnode = pe.Node( - # niu.IdentityInterface(fields=['qc', 'mosaic', 'out_group', 'out_dvars', 'out_fd']), - # name='outputnode', - # ) - sanitize = pe.Node( SanitizeImage( n_volumes_to_discard=0, @@ -141,14 +136,6 @@ def dmri_qc_workflow(name='dwiMRIQC'): mem_gb=mem_gb * 2.5, ) - # Shell-wise hmc not functional (yet?) - # hmc_shells = pe.MapNode( - # Volreg(args='-Fourier -twopass', zpad=4, outputtype='NIFTI_GZ'), - # name='hmc_shells', - # mem_gb=mem_gb * 2.5, - # iterfield=['in_file'], - # ) - # Calculate brainmask dmri_bmsk = dmri_bmsk_workflow(omp_nthreads=config.nipype.omp_nthreads) @@ -230,6 +217,7 @@ def dmri_qc_workflow(name='dwiMRIQC'): (dmri_bmsk, sp_mask, [('outputnode.out_mask', 'brain_mask')]), (dmri_bmsk, drift, [('outputnode.out_mask', 'brainmask_file')]), (drift, hmcwf, [('out_full_file', 'inputnode.in_file')]), + (meta, hmcwf, [('out_bvec_file', 'inputnode.in_bvec')]), (drift, averages, [('out_full_file', 'in_file')]), (drift, stddev, [('out_full_file', 'in_file')]), (shells, averages, [('b_masks', 'in_weights')]), @@ -251,7 +239,9 @@ def dmri_qc_workflow(name='dwiMRIQC'): (averages, iqms_wf, [(('out_file', _first), 'inputnode.in_b0')]), (sp_mask, iqms_wf, [('out_mask', 'inputnode.spikes_mask')]), (piesno, iqms_wf, [('sigma', 'inputnode.piesno_sigma')]), - (hmcwf, iqms_wf, [('outputnode.out_fd', 'inputnode.framewise_displacement')]), + (hmcwf, iqms_wf, [('outputnode.out_fd', 'inputnode.framewise_displacement'), + ('outputnode.out_bvec', 'inputnode.in_bvec_rotated'), + ('outputnode.out_bvec_diff', 'inputnode.in_bvec_diff')]), (dwimodel, iqms_wf, [('out_fa', 'inputnode.in_fa'), ('out_cfa', 'inputnode.in_cfa'), ('out_fa_nans', 'inputnode.in_fa_nans'), @@ -312,6 +302,8 @@ def compute_iqms(name='ComputeIQMs'): 'in_file', 'in_shells', 'in_bvec', + 'in_bvec_rotated', + 'in_bvec_diff', 'in_b0', 'in_fa', 'in_cfa', @@ -377,6 +369,8 @@ def compute_iqms(name='ComputeIQMs'): ('b_values', 'in_bval'), ('in_shells', 'in_shells'), ('in_bvec', 'in_bvec'), + ('in_bvec_rotated', 'in_bvec_rotated'), + ('in_bvec_diff', 'in_bvec_diff'), ('in_b0', 'in_b0'), ('brain_mask', 'brain_mask'), ('wm_mask', 'wm_mask'), @@ -426,12 +420,23 @@ def hmc_workflow(name='dMRI_HMC'): from nipype.algorithms.confounds import FramewiseDisplacement from nipype.interfaces.afni import Volreg + from mriqc.interfaces.diffusion import RotateVectors + mem_gb = config.workflow.biggest_file_gb workflow = pe.Workflow(name=name) - inputnode = pe.Node(niu.IdentityInterface(fields=['in_file', 'reference']), name='inputnode') - outputnode = pe.Node(niu.IdentityInterface(fields=['out_file', 'out_fd']), name='outputnode') + inputnode = pe.Node(niu.IdentityInterface(fields=[ + 'in_file', + 'reference', + 'in_bvec', + ]), name='inputnode') + outputnode = pe.Node(niu.IdentityInterface(fields=[ + 'out_file', + 'out_fd', + 'out_bvec', + 'out_bvec_diff', + ]), name='outputnode') # calculate hmc parameters hmc = pe.Node( @@ -440,6 +445,8 @@ def hmc_workflow(name='dMRI_HMC'): mem_gb=mem_gb * 2.5, ) + bvec_rot = pe.Node(RotateVectors(), name='bvec_rot') + # Compute the frame-wise displacement fdnode = pe.Node( FramewiseDisplacement( @@ -454,9 +461,14 @@ def hmc_workflow(name='dMRI_HMC'): workflow.connect([ (inputnode, hmc, [('in_file', 'in_file'), ('reference', 'basefile')]), + (inputnode, bvec_rot, [('in_bvec', 'in_file'), + ('reference', 'reference')]), (hmc, outputnode, [('out_file', 'out_file')]), (hmc, fdnode, [('oned_file', 'in_file')]), + (hmc, bvec_rot, [('oned_matrix_save', 'transforms')]), (fdnode, outputnode, [('out_file', 'out_fd')]), + (bvec_rot, outputnode, [('out_bvec', 'out_bvec'), + ('out_diff', 'out_bvec_diff')]), ]) # fmt: on return workflow