Merge pull request #200 from jhlegarreta/AddGradAngleComputationUtils

ENH: Add gradient encoding direction angle computation utils

src/eddymotion/math/tests/test_utils.py

@@ -21,8 +21,9 @@
 #     https://www.nipreps.org/community/licensing/
 #
 import numpy as np
+import pytest
 
-from eddymotion.math.utils import is_positive_definite
+from eddymotion.math.utils import compute_angle, is_positive_definite
 
 
 def test_is_positive_definite():
@@ -31,3 +32,17 @@ def test_is_positive_definite():
 
     matrix = np.array([[4, 1, 2], [1, -3, 1], [2, 1, 5]])
     assert not is_positive_definite(matrix)
+
+
+@pytest.mark.parametrize(
+    ("v1", "v2", "closest_polarity", "expected"),
+    [
+        ([1, 0, 0], [1, 0, 0], False, 0),
+        ([1, 0, 0], [0, 1, 0], False, np.pi / 2),
+        ([1, 0, 0], [-1 / np.sqrt(2), 0, 1 / np.sqrt(2)], False, np.pi * 3 / 4),
+        ([1, 0, 0], [-1 / np.sqrt(2), 0, 1 / np.sqrt(2)], True, np.pi / 4),
+    ],
+)
+def test_compute_angle(v1, v2, closest_polarity, expected):
+    obtained = compute_angle(v1, v2, closest_polarity=closest_polarity)
+    assert np.isclose(obtained, expected)

src/eddymotion/math/utils.py

@@ -29,7 +29,7 @@ def is_positive_definite(matrix):
     matrix and its conjugate transpose by performing the Cholesky decomposition.
     Parameters
     ----------
-    matrix : np.ndarray
+    matrix : :obj:`~numpy.ndarray`
         The matrix to check.
     Returns
     -------
@@ -43,3 +43,47 @@ def is_positive_definite(matrix):
     except np.linalg.LinAlgError:
         # Matrix is not positive definite
         return False
+
+
+def compute_angle(v1, v2, closest_polarity=False):
+    """Compute the angle between two vectors.
+
+    Parameters
+    ----------
+    v1 : :obj:`~numpy.ndarray`
+        First vector.
+    v2 : :obj:`~numpy.ndarray`
+        Second vector.
+    closest_polarity : :obj:`bool`
+        ``True`` to consider the smallest of the two angles between the crossing
+         lines resulting from reversing both vectors.
+
+    Returns
+    -------
+    :obj:`float`
+        The angle between the two vectors in radians.
+
+    Examples
+    --------
+    >>> compute_angle(
+    ...     np.array((1.0, 0.0, 0.0)),
+    ...     np.array((-1.0, 0.0, 0.0)),
+    ... )  # doctest: +ELLIPSIS
+    3.1415...
+    >>> compute_angle(
+    ...     np.array((1.0, 0.0, 0.0)),
+    ...     np.array((-1.0, 0.0, 0.0)),
+    ...     closest_polarity=True,
+    ... )
+    0.0
+
+    """
+
+    cosine_angle = (v1 / np.linalg.norm(v1)) @ (v2 / np.linalg.norm(v2))
+    # Clip values to handle numerical errors
+    cosine_angle = np.clip(
+        np.abs(cosine_angle) if closest_polarity else cosine_angle,
+        -1.0,
+        1.0,
+    )
+    return np.arccos(cosine_angle)

src/eddymotion/model/gradient_utils.py

@@ -0,0 +1,75 @@
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+#
+# Copyright 2024 The NiPreps Developers <[email protected]>
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+# We support and encourage derived works from this project, please read
+# about our expectations at
+#
+#     https://www.nipreps.org/community/licensing/
+#
+import numpy as np
+
+
+def compute_pairwise_angles(bvecs, closest_polarity):
+    r"""Compute pairwise angles across diffusion gradient encoding directions.
+
+    Following [Andersson15]_, it computes the smallest of the angles between
+    each pair if ``closest_polarity`` is ``True``, i.e.
+
+    .. math::
+
+        \theta(\mathbf{g}, \mathbf{g'}) = \arccos(\abs{\langle \mathbf{g}, \mathbf{g'} \rangle})
+
+    Parameters
+    ----------
+    bvecs : :obj:`~numpy.ndarray`
+        Diffusion gradient encoding directions in FSL format.
+    closest_polarity : :obj:`bool`
+        ``True`` to consider the smallest of the two angles between the crossing
+         lines resulting from reversing each vector pair.
+
+    Returns
+    -------
+    :obj:`~numpy.ndarray`
+        Pairwise angles across diffusion gradient encoding directions.
+
+    Examples
+    --------
+    >>> compute_pairwise_angles(
+    ...     ((1.0, -1.0), (0.0, 0.0), (0.0, 0.0)),
+    ...     False,
+    ... )[0, 1]  # doctest: +ELLIPSIS
+    3.1415...
+    >>> compute_pairwise_angles(
+    ...     ((1.0, -1.0), (0.0, 0.0), (0.0, 0.0)),
+    ...     True,
+    ... )[0, 1]
+    0.0
+
+    References
+    ----------
+    .. [Andersson15] J. L. R. Andersson. et al., An integrated approach to
+       correction for off-resonance effects and subject movement in diffusion MR
+       imaging, NeuroImage 125 (2016) 1063–1078
+    """
+
+    if np.shape(bvecs)[0] != 3:
+        raise ValueError(f"bvecs must be of shape (3, N). Found: {bvecs.shape}")
+
+    # Ensure b-vectors are unit-norm
+    bvecs = np.array(bvecs) / np.linalg.norm(bvecs, axis=0)
+    cosines = np.clip(bvecs.T @ bvecs, -1.0, 1.0)
+    return np.arccos(np.abs(cosines) if closest_polarity else cosines)

src/eddymotion/model/tests/test_gradient_utils.py

@@ -0,0 +1,62 @@
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+#
+# Copyright 2024 The NiPreps Developers <[email protected]>
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+# We support and encourage derived works from this project, please read
+# about our expectations at
+#
+#     https://www.nipreps.org/community/licensing/
+#
+import numpy as np
+
+from eddymotion.model.gradient_utils import compute_pairwise_angles
+
+
+def test_compute_pairwise_angles():
+    # No need to use normalized vectors: compute_angle takes care of dealing
+    # with it.
+    # The last vector serves as a case where e.g. the angle between the first
+    # vector and the last one is 135, and the method yielding the smallest
+    # resulting angle between the crossing lines (45 vs 135)
+    bvecs = np.array(
+        [
+            [1, 0, 0, 1, 1, 0, -1],
+            [0, 1, 0, 1, 0, 1, 0],
+            [0, 0, 1, 0, 1, 1, 1],
+        ]
+    )
+
+    expected = np.array(
+        [
+            [0.0, np.pi / 2, np.pi / 2, np.pi / 4, np.pi / 4, np.pi / 2, np.pi / 4],
+            [np.pi / 2, 0.0, np.pi / 2, np.pi / 4, np.pi / 2, np.pi / 4, np.pi / 2],
+            [np.pi / 2, np.pi / 2, 0.0, np.pi / 2, np.pi / 4, np.pi / 4, np.pi / 4],
+            [np.pi / 4, np.pi / 4, np.pi / 2, 0.0, np.pi / 3, np.pi / 3, np.pi / 3],
+            [np.pi / 4, np.pi / 2, np.pi / 4, np.pi / 3, 0.0, np.pi / 3, np.pi / 2],
+            [np.pi / 2, np.pi / 4, np.pi / 4, np.pi / 3, np.pi / 3, 0.0, np.pi / 3],
+            [np.pi / 4, np.pi / 2, np.pi / 4, np.pi / 3, np.pi / 2, np.pi / 3, 0.0],
+        ]
+    )
+
+    smallest = True
+    obtained = compute_pairwise_angles(bvecs, smallest)
+
+    # Expect N*N elements
+    assert bvecs.shape[-1] ** 2 == np.prod(obtained.shape)
+    assert obtained.shape == expected.shape
+    # Check that the matrix is symmetric
+    assert np.allclose(expected, expected.T)
+    np.testing.assert_array_almost_equal(obtained, expected, decimal=2)