Merge pull request #903 from EmmaRenauld/test_gradients_utils

Add gradient_utils unit tests

scilpy/gradients/tests/test_gradients_utils.py

@@ -1,10 +1,44 @@
 # -*- coding: utf-8 -*-
+import numpy as np
+
+from scilpy.gradients.bvec_bval_tools import is_normalized_bvecs
+from scilpy.gradients.utils import random_uniform_on_sphere, \
+    get_new_gtab_order
 
 
 def test_random_uniform_on_sphere():
-    pass
+    bvecs = random_uniform_on_sphere(10)
+
+    # Confirm that they are unit vectors.
+    assert is_normalized_bvecs(bvecs)
+
+    # Can't check much more than that. Supervising that no co-linear vectors.
+    # (Each pair of vector should have at least some angle in-between)
+    # We also tried to check if the closest neighbor to each vector is more or
+    # less always with the same angle, but it is very variable.
+    min_expected_angle = 1.0
+    smallests = []
+    for i in range(10):
+        angles = np.rad2deg(np.arccos(np.dot(bvecs[i, :], bvecs.T)))
+        # Smallest, except 0 (with itself). Sometimes this is nan.
+        smallests.append(np.nanmin(angles[angles > 1e-5]))
+    assert np.all(np.asarray(smallests) > min_expected_angle)
+
+
+def test_get_new_gtab_order():
+    # Using N=4 vectors
+    philips_table = np.asarray([[1, 1, 1, 1],
+                                [2, 2, 2, 2],
+                                [3, 3, 3, 3],
+                                [4, 4, 4, 4]])
+    dwi = np.ones((10, 10, 10, 4))
+    bvecs = np.asarray([[3, 3, 3],
+                        [4, 4, 4],
+                        [2, 2, 2],
+                        [1, 1, 1]])
+    bvals = np.asarray([3, 4, 2, 1])
 
+    order = get_new_gtab_order(philips_table, dwi, bvals, bvecs)
 
-def test_get_new_order_philips():
-    # Needs dwi files - philips version
-    pass
+    assert np.array_equal(bvecs[order, :], philips_table[:, 0:3])
+    assert np.array_equal(bvals[order], philips_table[:, 3])

scilpy/gradients/utils.py

@@ -6,7 +6,11 @@
 def random_uniform_on_sphere(nb_vectors):
     """
     Creates a set of K pseudo-random unit vectors, following a uniform
-    distribution on the sphere.
+    distribution on the sphere. Reference: Emmanuel Caruyer's code
+    (https://github.com/ecaruyer).
+
+    This is not intended to create a perfect result. It's usually the
+    initialization step of a repulsion strategy.
 
     Parameters
     ----------
@@ -15,14 +19,21 @@ def random_uniform_on_sphere(nb_vectors):
 
     Returns
     -------
-    bvecs: nd.array
-        pseudo-random unit vector
+    bvecs: nd.array of shape (nb_vectors, 3)
+        Pseudo-random unit vectors
     """
+    # Note. Caruyer's docstring says it's a uniform on the half-sphere, but
+    # plotted a few results: it is one the whole sphere.
     phi = 2 * np.pi * np.random.rand(nb_vectors)
 
     r = 2 * np.sqrt(np.random.rand(nb_vectors))
     theta = 2 * np.arcsin(r / 2)
 
+    # See here:
+    # https://www.bogotobogo.com/Algorithms/uniform_distribution_sphere.php
+    # They seem to be using something like this instead:
+    # theta = np.arccos(2 * np.random.rand(nb_vectors) - 1.0)
+
     bvecs = np.zeros((nb_vectors, 3))
     bvecs[:, 0] = np.sin(theta) * np.cos(phi)
     bvecs[:, 1] = np.sin(theta) * np.sin(phi)
@@ -31,54 +42,59 @@ def random_uniform_on_sphere(nb_vectors):
     return bvecs
 
 
-def get_new_order_philips(philips_table, dwi, bvals, bvecs):
+def get_new_gtab_order(ref_gradient_table, dwi, bvals, bvecs):
     """
-    Reorder bval and bvec files based on the philips table.
+    Find the sorting order that could be applied to the bval and bvec files to
+    obtain the same order as in the reference gradient table.
+
+    This is mostly useful to reorder bval and bvec files in the order they were
+    acquired by the Philips scanner (before version 5.6).
 
     Parameters
     ----------
-    philips_table: nd.array
-        Philips gradient table
-    dwis: nibabel
-        dwis
+    ref_gradient_table: nd.array
+        Gradient table, of shape (N, 4). It will use as reference for the
+        ordering of b-vectors.
+        Ex: Could be the result of scil_gradients_generate_sampling.py
+    dwi: nibabel image
+        dwi of shape (x, y, z, N). Only used to confirm the dwi's shape.
     bvals : array, (N,)
-        bvals
+        bvals that need to be reordered.
     bvecs : array, (N, 3)
-        bvecs
+        bvecs that need to be reordered.
 
     Returns
     -------
     new_index: nd.array
         New index to reorder bvals/bvec
     """
-    # Check number of gradients, bvecs, bvals, dwi and oTable
-    if len(bvecs) != dwi.shape[3] or len(bvals) != len(philips_table):
-        raise ValueError('bvec/bval/dwi and original table \
-            does not contain the same number of gradients')
-
-    # Check bvals
-    philips_bval = np.unique(philips_table[:, 3])
+    if not (len(bvecs) == dwi.shape[3] == len(bvals) ==
+            len(ref_gradient_table)):
+        raise ValueError('bvec/bval/dwi and reference table do not contain '
+                         'the same number of gradients.')
 
-    philips_dwi_shells = philips_bval[philips_bval > 1]
-    philips_b0s = philips_bval[philips_bval < 1]
+    ref_bval = np.unique(ref_gradient_table[:, 3])
+    ref_dwi_shells = ref_bval[ref_bval > 1]
+    ref_b0s = ref_bval[ref_bval < 1]
 
     dwi_shells = np.unique(bvals[bvals > 1])
     b0s = np.unique(bvals[bvals < 1])
 
-    if len(philips_dwi_shells) != len(dwi_shells) or\
-       len(philips_b0s) != len(b0s):
-        raise ValueError('bvec/bval/dwi and original table\
-                          does not contain the same shells')
+    if len(ref_dwi_shells) != len(dwi_shells) or \
+       len(ref_b0s) != len(b0s):
+        raise ValueError('bvec/bval/dwi and reference table do not contain '
+                         'the same shells.')
 
     new_index = np.zeros(bvals.shape)
 
-    for nbval in philips_bval:
+    for nbval in ref_bval:
         curr_bval = np.where(bvals == nbval)[0]
-        curr_bval_table = np.where(philips_table[:, 3] == nbval)[0]
+        curr_bval_table = np.where(ref_gradient_table[:, 3] == nbval)[0]
 
         if len(curr_bval) != len(curr_bval_table):
-            raise ValueError('bval/bvec and orginal table does not contain \
-                the same number of gradients for shell {0}'.format(curr_bval))
+            raise ValueError('bval/bvec and orginal table do not contain '
+                             'the same number of gradients for shell {}.'
+                             .format(curr_bval))
 
         new_index[curr_bval_table] = curr_bval
 

scripts/scil_dwi_reorder_philips.py

@@ -18,7 +18,7 @@
 import nibabel as nib
 import numpy as np
 
-from scilpy.gradients.utils import get_new_order_philips
+from scilpy.gradients.utils import get_new_gtab_order
 from scilpy.io.utils import (add_overwrite_arg,
                              add_verbose_arg,
                              assert_inputs_exist,
@@ -88,7 +88,7 @@ def main():
     bvals, bvecs = read_bvals_bvecs(args.in_bval, args.in_bvec)
     dwi = nib.load(args.in_dwi)
 
-    new_index = get_new_order_philips(philips_table, dwi, bvals, bvecs)
+    new_index = get_new_gtab_order(philips_table, dwi, bvals, bvecs)
     bvecs = bvecs[new_index]
     bvals = bvals[new_index]