diff --git a/dmriprep/interfaces/vectors.py b/dmriprep/interfaces/vectors.py index 6d9da11b..223fc91f 100644 --- a/dmriprep/interfaces/vectors.py +++ b/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 diff --git a/dmriprep/utils/vectors.py b/dmriprep/utils/vectors.py index d9284ef3..1dd2faea 100644 --- a/dmriprep/utils/vectors.py +++ b/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]))