[ENH] Reorient bvecs from affs

dmriprep/interfaces/vectors.py

@@ -6,7 +6,8 @@
     SimpleInterface, BaseInterfaceInputSpec, TraitedSpec,
     File, traits, isdefined
 )
-from ..utils.vectors import DiffusionGradientTable, B0_THRESHOLD, BVEC_NORM_EPSILON
+from ..utils.vectors import DiffusionGradientTable, reorient_vecs_from_ras_b, \
+    B0_THRESHOLD, BVEC_NORM_EPSILON
 
 
 class _CheckGradientTableInputSpec(BaseInterfaceInputSpec):
@@ -97,3 +98,62 @@ def _undefined(objekt, name, default=None):
     if not isdefined(value):
         return default
     return value
+
+
+class _ReorientVectorsInputSpec(BaseInterfaceInputSpec):
+    in_rasb = File(exists=True)
+    affines = traits.Array()
+    b0_threshold = traits.Float(B0_THRESHOLD, usedefault=True)
+
+
+class _ReorientVectorsOutputSpec(TraitedSpec):
+    out_rasb = File(exists=True)
+
+
+class ReorientVectors(SimpleInterface):
+    """
+    Reorient Vectors
+
+    Example
+    -------
+
+    >>> os.chdir(tmpdir)
+    >>> oldrasb = str(data_dir / 'dwi.tsv')
+    >>> oldrasb_mat = np.loadtxt(str(data_dir / 'dwi.tsv'), skiprows=1)
+    >>> # The simple case: all affines are identity
+    >>> affine_list = np.zeros((len(oldrasb_mat[:, 3][oldrasb_mat[:, 3] != 0]), 4, 4))
+    >>> for i in range(4):
+    ...     affine_list[:, i, i] = 1
+    ...     reor_vecs = ReorientVectors()
+    >>> reor_vecs = ReorientVectors()
+    >>> reor_vecs.inputs.affines = affine_list
+    >>> reor_vecs.inputs.in_rasb = oldrasb
+    >>> res = reor_vecs.run()
+    >>> out_rasb = res.outputs.out_rasb
+    >>> out_rasb_mat = np.loadtxt(out_rasb, skiprows=1)
+    >>> npt.assert_equal(oldrasb_mat, out_rasb_mat)
+    True
+    """
+
+    input_spec = _ReorientVectorsInputSpec
+    output_spec = _ReorientVectorsOutputSpec
+
+    def _run_interface(self, runtime):
+        rasb_file = _undefined(self.inputs, 'in_rasb')
+
+        reor_table = reorient_vecs_from_ras_b(
+            rasb_file=rasb_file,
+            affines=self.inputs.affines,
+            b0_threshold=self.inputs.b0_threshold,
+        )
+
+        cwd = Path(runtime.cwd).absolute()
+        reor_rasb_file = fname_presuffix(
+            self.inputs.in_rasb, use_ext=False, suffix='_reoriented.tsv',
+            newpath=str(cwd))
+        np.savetxt(str(reor_rasb_file), reor_table,
+                   delimiter='\t', header='\t'.join('RASB'),
+                   fmt=['%.8f'] * 3 + ['%g'])
+
+        self._results['out_rasb'] = reor_rasb_file
+        return runtime

dmriprep/utils/vectors.py

@@ -375,3 +375,50 @@ def bvecs2ras(affine, bvecs, norm=True, bvec_norm_epsilon=0.2):
         rotated_bvecs[~b0s] /= norms_bvecs[~b0s, np.newaxis]
         rotated_bvecs[b0s] = np.zeros(3)
     return rotated_bvecs
+
+
+def reorient_vecs_from_ras_b(rasb_file, affines, b0_threshold=B0_THRESHOLD):
+    """
+    Reorient the vectors from a rasb .tsv file.
+    When correcting for motion, rotation of the diffusion-weighted volumes
+    might cause systematic bias in rotationally invariant measures, such as FA
+    and MD, and also cause characteristic biases in tractography, unless the
+    gradient directions are appropriately reoriented to compensate for this
+    effect [Leemans2009]_.
+
+    Parameters
+    ----------
+    rasb_file : str or os.pathlike
+        File path to a RAS-B gradient table.
+
+    affines : list or ndarray of shape (n, 4, 4) or (n, 3, 3)
+        Each entry in this list or array contain either an affine
+        transformation (4,4) or a rotation matrix (3, 3).
+        In both cases, the transformations encode the rotation that was applied
+        to the image corresponding to one of the non-zero gradient directions.
+
+    Returns
+    -------
+    Gradients : ndarray of shape (4, n)
+        A reoriented ndarray where the first three columns correspond to each of
+        x, y, z directions of the bvecs, in R-A-S image orientation.
+
+    """
+    from dipy.core.gradients import gradient_table_from_bvals_bvecs, reorient_bvecs
+
+    ras_b_mat = np.genfromtxt(rasb_file, delimiter='\t')
+
+    # Verify that number of non-B0 volumes corresponds to the number of affines.
+    # If not, raise an error.
+    if len(ras_b_mat[:, 3][ras_b_mat[:, 3] != 0]) != len(affines):
+        raise ValueError('Number of affine transformations must match number of '
+                         'non-zero gradients')
+
+    # Build gradient table object
+    gt = gradient_table_from_bvals_bvecs(ras_b_mat[:, 3], ras_b_mat[:, 0:3],
+                                         b0_threshold=b0_threshold)
+
+    # Reorient table
+    new_gt = reorient_bvecs(gt, affines)
+
+    return np.hstack((new_gt.bvecs, new_gt.bvals[..., np.newaxis]))