diff --git a/scilpy/gradients/bvec_bval_tools.py b/scilpy/gradients/bvec_bval_tools.py
index 2bcf53848..fbc172395 100644
--- a/scilpy/gradients/bvec_bval_tools.py
+++ b/scilpy/gradients/bvec_bval_tools.py
@@ -4,6 +4,9 @@
 from enum import Enum
 
 from dipy.core.gradients import get_bval_indices
+from dipy.core.sphere import Sphere
+from dipy.data import get_sphere
+from scipy.spatial import KDTree
 import numpy as np
 
 DEFAULT_B0_THRESHOLD = 20
@@ -313,3 +316,36 @@ def round_bvals_to_shell(bvals, shells_to_extract, tol=20):
                          '''.format(bvals[modified.astype(bool)]))
 
     return new_bvals
+
+
+def compare_bvecs_to_sphere(bvecs):
+    """
+    Finds the maximum angle between any sphere vertice and its nearest bvec.
+    """
+    tree = KDTree(normalize_bvecs(bvecs))
+
+    sphere = get_sphere()
+
+    print(bvecs[:10])
+    print(sphere.vertices[:10])
+
+    distances, ids = tree.query(sphere.vertices,
+                              workers=-1)
+    print(np.ones((len(sphere.vertices))))
+    print(distances, ids)
+
+    most_lonely_vertice_id = np.argmax(distances)
+    most_lonely_vertice = sphere.vertices[most_lonely_vertice_id]
+    max_nearest_bvec_id = ids[most_lonely_vertice_id]
+    max_nearest_bvec = bvecs[max_nearest_bvec_id]
+
+    print(most_lonely_vertice)
+    print(np.max(distances))
+    max_angle = np.arccos(
+        np.clip([np.dot(max_nearest_bvec, most_lonely_vertice)], [-1], [1]))
+
+    print(max_angle)
+
+    # TODO: Add energy score
+
+    return np.rad2deg(max_angle)
diff --git a/scripts/scil_gradients_validate_sphere.py b/scripts/scil_gradients_validate_sphere.py
new file mode 100755
index 000000000..84071d99f
--- /dev/null
+++ b/scripts/scil_gradients_validate_sphere.py
@@ -0,0 +1,51 @@
+#! /usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Compares a gradient table with the Dipy sphere and outputs the maximal angle
+between any sphere vertice and its closest bvec.
+"""
+
+import argparse
+import logging
+
+from dipy.io.gradients import read_bvals_bvecs
+
+from scilpy.io.utils import (add_overwrite_arg, assert_inputs_exist,
+                             add_json_args, add_verbose_arg,
+                             assert_outputs_exist)
+from scilpy.gradients.bvec_bval_tools import compare_bvecs_to_sphere
+
+
+def _build_arg_parser():
+    p = argparse.ArgumentParser(formatter_class=argparse.RawTextHelpFormatter,
+                                description=__doc__)
+
+    p.add_argument('in_bvec',
+                   help='Path to bvec file.')
+    p.add_argument('out_metrics',
+                   help='Path to a json file in which the validation results\n'
+                        'will be saved.')
+
+    add_json_args(p)
+    add_verbose_arg(p)
+    add_overwrite_arg(p)
+    return p
+
+
+def main():
+    parser = _build_arg_parser()
+    args = parser.parse_args()
+    logging.getLogger().setLevel(logging.getLevelName(args.verbose))
+
+    assert_inputs_exist(parser, [args.in_bvec])
+    assert_outputs_exist(parser, args, args.out_metrics)
+
+    logging.debug("Loading bvecs")
+    _, bvecs = read_bvals_bvecs(None, args.in_bvec)
+
+    max_angle = compare_bvecs_to_sphere(bvecs)
+
+    print(max_angle)
+
+if __name__ == "__main__":
+    main()