Merge pull request #669 from arokem/reorient-bvecs

Function to reorient gradient directions according to moco parameters

dipy/core/gradients.py

@@ -3,6 +3,11 @@
 from ..utils.six import string_types
 
 import numpy as np
+try:
+    from scipy.linalg import polar
+except ImportError:   # Some elderly scipy doesn't have polar
+    from dipy.fixes.scipy import polar
+from scipy.linalg import inv
 
 from ..io import gradients as io
 from .onetime import auto_attr
@@ -88,7 +93,7 @@ def info(self):
 
 
 def gradient_table_from_bvals_bvecs(bvals, bvecs, b0_threshold=0, atol=1e-2,
-                                  **kwargs):
+                                    **kwargs):
     """Creates a GradientTable from a bvals array and a bvecs array
 
     Parameters
@@ -220,9 +225,9 @@ def gradient_table(bvals, bvecs=None, big_delta=None, small_delta=None,
     # If you provided strings with full paths, we go and load those from
     # the files:
     if isinstance(bvals, string_types):
-          bvals, _ = io.read_bvals_bvecs(bvals, None)
+        bvals, _ = io.read_bvals_bvecs(bvals, None)
     if isinstance(bvecs, string_types):
-          _, bvecs = io.read_bvals_bvecs(None, bvecs)
+        _, bvecs = io.read_bvals_bvecs(None, bvecs)
 
     bvals = np.asarray(bvals)
     # If bvecs is None we expect bvals to be an (N, 4) or (4, N) array.
@@ -238,9 +243,66 @@ def gradient_table(bvals, bvecs=None, big_delta=None, small_delta=None,
                              " array containing both bvals and bvecs")
     else:
         bvecs = np.asarray(bvecs)
-        if (bvecs.shape[1] > bvecs.shape[0])  and bvecs.shape[0]>1:
+        if (bvecs.shape[1] > bvecs.shape[0]) and bvecs.shape[0] > 1:
             bvecs = bvecs.T
     return gradient_table_from_bvals_bvecs(bvals, bvecs, big_delta=big_delta,
                                            small_delta=small_delta,
                                            b0_threshold=b0_threshold,
                                            atol=atol)
+
+
+def reorient_bvecs(gtab, affines):
+    """Reorient the directions in a GradientTable.
+
+    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
+    ----------
+    gtab : GradientTable
+        The nominal gradient table with which the data were acquired.
+    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
+        (ordered according to their order in `gtab.bvecs[~gtab.b0s_mask]`)
+
+    Returns
+    -------
+    gtab : a GradientTable class instance with the reoriented directions
+
+    References
+    ----------
+    .. [Leemans2009] The B-Matrix Must Be Rotated When Correcting for
+       Subject Motion in DTI Data. Leemans, A. and Jones, D.K. (2009).
+       MRM, 61: 1336-1349
+    """
+    new_bvecs = gtab.bvecs[~gtab.b0s_mask]
+
+    if new_bvecs.shape[0] != len(affines):
+        e_s = "Number of affine transformations must match number of "
+        e_s += "non-zero gradients"
+        raise ValueError(e_s)
+
+    for i, aff in enumerate(affines):
+        if aff.shape == (4, 4):
+            # This must be an affine!
+            # Remove the translation component:
+            aff_no_trans = aff[:3, :3]
+            # Decompose into rotation and scaling components:
+            R, S = polar(aff_no_trans)
+        elif aff.shape == (3, 3):
+            # We assume this is a rotation matrix:
+            R = aff
+        Rinv = inv(R)
+        # Apply the inverse of the rotation to the corresponding gradient
+        # direction:
+        new_bvecs[i] = np.dot(Rinv, new_bvecs[i])
+
+    return_bvecs = np.zeros(gtab.bvecs.shape)
+    return_bvecs[~gtab.b0s_mask] = new_bvecs
+    return gradient_table(gtab.bvals, return_bvecs)

dipy/core/tests/test_gradients.py

@@ -1,11 +1,12 @@
-from nose.tools import assert_true
+from nose.tools import assert_true, assert_raises
 
 import numpy as np
 import numpy.testing as npt
 
 from dipy.data import get_data
 from dipy.core.gradients import (gradient_table, GradientTable,
-                                 gradient_table_from_bvals_bvecs)
+                                 gradient_table_from_bvals_bvecs,
+                                 reorient_bvecs)
 from dipy.io.gradients import read_bvals_bvecs
 
 
@@ -39,6 +40,8 @@ def test_btable_prepare():
     bt4 = gradient_table(btab.T)
     npt.assert_array_equal(bt4.bvecs, bvecs)
     npt.assert_array_equal(bt4.bvals, bvals)
+    # Test for proper inputs (expects either bvals/bvecs or 4 by n):
+    assert_raises(ValueError, gradient_table, bvecs)
 
 
 def test_GradientTable():
@@ -183,6 +186,69 @@ def test_qvalues():
     bt = gradient_table(bvals, bvecs, big_delta=8, small_delta=6)
     npt.assert_almost_equal(bt.qvals, qvals)
 
+
+def test_reorient_bvecs():
+    sq2 = np.sqrt(2) / 2
+    bvals = np.concatenate([[0], np.ones(6) * 1000])
+    bvecs = np.array([[0, 0, 0],
+                      [1, 0, 0],
+                      [0, 1, 0],
+                      [0, 0, 1],
+                      [sq2, sq2, 0],
+                      [sq2, 0, sq2],
+                      [0, sq2, sq2]])
+
+    gt = gradient_table_from_bvals_bvecs(bvals, bvecs, b0_threshold=0)
+    # The simple case: all affines are identity
+    affs = np.zeros((6, 4, 4))
+    for i in range(4):
+        affs[:, i, i] = 1
+
+    # We should get back the same b-vectors
+    new_gt = reorient_bvecs(gt, affs)
+    npt.assert_equal(gt.bvecs, new_gt.bvecs)
+
+    # Now apply some rotations
+    rotation_affines = []
+    rotated_bvecs = bvecs[:]
+    for i in np.where(~gt.b0s_mask)[0]:
+        rot_ang = np.random.rand()
+        cos_rot = np.cos(rot_ang)
+        sin_rot = np.sin(rot_ang)
+        rotation_affines.append(np.array([[1, 0, 0, 0],
+                                         [0, cos_rot, -sin_rot, 0],
+                                         [0, sin_rot, cos_rot, 0],
+                                         [0, 0, 0, 1]]))
+        rotated_bvecs[i] = np.dot(rotation_affines[-1][:3, :3],
+                                  bvecs[i])
+
+    # Copy over the rotation affines
+    full_affines = rotation_affines[:]
+    # And add some scaling and translation to each one to make this harder
+    for i in range(len(full_affines)):
+        full_affines[i] = np.dot(full_affines[i],
+                                 np.array([[2.5, 0, 0, -10],
+                                           [0, 2.2, 0, 20],
+                                           [0, 0, 1, 0],
+                                           [0, 0, 0, 1]]))
+
+    gt_rot = gradient_table_from_bvals_bvecs(bvals,
+                                             rotated_bvecs, b0_threshold=0)
+    new_gt = reorient_bvecs(gt_rot, full_affines)
+    # At the end of all this, we should be able to recover the original
+    # vectors
+    npt.assert_almost_equal(gt.bvecs, new_gt.bvecs)
+
+    # We should be able to pass just the 3-by-3 rotation components to the same
+    # effect
+    new_gt = reorient_bvecs(gt_rot, np.array(rotation_affines)[:, :3, :3])
+    npt.assert_almost_equal(gt.bvecs, new_gt.bvecs)
+
+    # Verify that giving the wrong number of affines raises an error:
+    full_affines.append(np.zeros((4, 4)))
+    assert_raises(ValueError, reorient_bvecs, gt_rot, full_affines)
+
+
 if __name__ == "__main__":
     from numpy.testing import run_module_suite
     run_module_suite()

dipy/fixes/scipy.py

@@ -0,0 +1,95 @@
+from __future__ import division, print_function, absolute_import
+
+import numpy as np
+from scipy.linalg import svd
+
+
+__all__ = ['polar']
+
+
+def polar(a, side="right"):
+    """
+    Compute the polar decomposition.
+
+    Returns the factors of the polar decomposition [1]_ `u` and `p` such
+    that ``a = up`` (if `side` is "right") or ``a = pu`` (if `side` is
+    "left"), where `p` is positive semidefinite.  Depending on the shape
+    of `a`, either the rows or columns of `u` are orthonormal.  When `a`
+    is a square array, `u` is a square unitary array.  When `a` is not
+    square, the "canonical polar decomposition" [2]_ is computed.
+
+    Parameters
+    ----------
+    a : array_like, shape (m, n).
+        The array to be factored.
+    side : {'left', 'right'}, optional
+        Determines whether a right or left polar decomposition is computed.
+        If `side` is "right", then ``a = up``.  If `side` is "left",  then
+        ``a = pu``.  The default is "right".
+
+    Returns
+    -------
+    u : ndarray, shape (m, n)
+        If `a` is square, then `u` is unitary.  If m > n, then the columns
+        of `a` are orthonormal, and if m < n, then the rows of `u` are
+        orthonormal.
+    p : ndarray
+        `p` is Hermitian positive semidefinite.  If `a` is nonsingular, `p`
+        is positive definite.  The shape of `p` is (n, n) or (m, m), depending
+        on whether `side` is "right" or "left", respectively.
+
+    References
+    ----------
+    .. [1] R. A. Horn and C. R. Johnson, "Matrix Analysis", Cambridge
+           University Press, 1985.
+    .. [2] N. J. Higham, "Functions of Matrices: Theory and Computation",
+           SIAM, 2008.
+
+    Notes
+    -----
+    Copyright (c) 2001, 2002 Enthought, Inc.
+    All rights reserved.
+
+    Copyright (c) 2003-2012 SciPy Developers.
+    All rights reserved.
+
+    Redistribution and use in source and binary forms, with or without
+    modification, are permitted provided that the following conditions are met:
+
+    a. Redistributions of source code must retain the above copyright notice,
+     this list of conditions and the following disclaimer.
+    b. Redistributions in binary form must reproduce the above copyright
+     notice, this list of conditions and the following disclaimer in the
+     documentation and/or other materials provided with the distribution.
+    c. Neither the name of Enthought nor the names of the SciPy Developers
+     may be used to endorse or promote products derived from this software
+     without specific prior written permission.
+
+
+    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+    AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+    IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+    ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS
+    BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
+    OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
+    THE POSSIBILITY OF SUCH DAMAGE.
+    """
+    if side not in ['right', 'left']:
+        raise ValueError("`side` must be either 'right' or 'left'")
+    a = np.asarray(a)
+    if a.ndim != 2:
+        raise ValueError("`a` must be a 2-D array.")
+
+    w, s, vh = svd(a, full_matrices=False)
+    u = w.dot(vh)
+    if side == 'right':
+        # a = up
+        p = (vh.T.conj() * s).dot(vh)
+    else:
+        # a = pu
+        p = (w * s).dot(w.T.conj())
+    return u, p