enh: add b-vector angular deviations as IQMs

References: #1216.

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.

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