NF - Added cython utility functions

dipy/utils/fast_numpy.pxd

@@ -4,6 +4,13 @@
 
 cimport numpy as np
 
+from libc.stdlib cimport rand
+from libc.math cimport sqrt
+
+
+cdef extern from "limits.h":
+    int INT_MAX
+
 # Replaces a numpy.searchsorted(arr, number, 'right')
 cdef int where_to_insert(
         np.float_t* arr,
@@ -22,3 +29,21 @@ cdef void copy_point(
 cdef void scalar_muliplication_point(
         double * a,
         double scalar) nogil
+
+cpdef double random() nogil
+
+cpdef double norm(
+        double[:] v) nogil
+
+cpdef double dot(
+        double[:] v1, 
+        double[:] v2) nogil
+
+cpdef void normalize(
+        double[:] v) nogil
+
+cpdef void cross(
+        double[:] out, 
+        double[:] v1, 
+        double[:] v2) nogil
+

dipy/utils/fast_numpy.pyx

@@ -45,3 +45,90 @@ cdef void scalar_muliplication_point(
         int i = 0
     for i in range(3):
         a[i] *= scalar
+
+
+cpdef double random() nogil:
+    """Sample a random number between 0 and 1.
+
+    Returns
+    -------
+    _ : double
+        random number
+    """
+    return rand() / float(INT_MAX)
+
+
+cpdef double norm(double[:] v) nogil:
+    """Compute the vector norm.
+
+    Parameters
+    ----------
+    v : double[3]
+        input vector.
+
+    Returns
+    -------
+    _ : double
+        norm of the input vector.
+    """
+    return sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2])
+
+
+cpdef double dot(double[:] v1, double[:] v2) nogil:
+    """Compute vectors dot product.
+
+    Parameters
+    ----------
+    v1 : double[3]
+        input vector 1.
+    v2 : double[3]
+        input vector 2.
+
+    Returns
+    -------
+    _ : double
+        dot product of input vectors.
+    """
+    return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]
+
+
+cpdef void normalize(double[:] v) nogil:
+    """Normalize the vector.
+
+    Parameters
+    ----------
+    v : double[3]
+        input vector
+
+    Notes
+    -----
+    Overwrites the first argument.
+
+    """
+    cdef double scale = 1.0 / norm(v)
+    v[0] = v[0] * scale
+    v[1] = v[1] * scale
+    v[2] = v[2] * scale
+
+
+cpdef void cross(double[:] out, double[:] v1, double[:] v2) nogil:
+    """Compute vectors cross product. 
+
+    Parameters
+    ----------
+    out : double[3]
+        output vector.
+    v1 : double[3]
+        input vector 1.
+    v2 : double[3]
+        input vector 2.
+
+    Notes
+    -----
+    Overwrites the first argument.
+
+    """
+    out[0] = v1[1] * v2[2] - v1[2] * v2[1]
+    out[1] = v1[2] * v2[0] - v1[0] * v2[2]
+    out[2] = v1[0] * v2[1] - v1[1] * v2[0]
+

dipy/utils/tests/test_fast_numpy.py

@@ -0,0 +1,72 @@
+import timeit
+
+import numpy as np
+from numpy.testing import assert_, assert_almost_equal, assert_array_equal 
+from dipy.utils.fast_numpy import random, norm, normalize, dot, cross
+
+
+def test_random():
+    # Test that random numbers are between 0 and 1.
+    for _ in range(10):
+        vec = random()
+        assert_(vec <= 1 and vec >= 0)
+
+
+def test_norm():
+    # Test that the norm equal numpy norm.
+    for _ in range(10):
+        vec = np.random.random(3)
+        assert_almost_equal(norm(vec), np.linalg.norm(vec))
+
+
+def test_normalize():
+    # Test that the normalize vector as a norm of 1.
+    for _ in range(10):
+        vec = np.random.random(3)
+        normalize(vec)
+        assert_almost_equal(np.linalg.norm(vec), 1)
+
+def test_dot():
+    # Test that dot is faster and equal to numpy.dot.
+    for _ in range(10):
+        vec1 = np.random.random(3)
+        vec2 = np.random.random(3)
+        assert_almost_equal(dot(vec1, vec2), np.dot(vec1, vec2))
+
+    vec1 = np.random.random(3)
+    vec2 = np.random.random(3)
+
+    def __dot():
+        dot(vec1, vec2)
+
+    def __npdot():
+        np.dot(vec1, vec2)
+
+    number=100000
+    time_dot = timeit.timeit(__dot, number=number)
+    time_npdot = timeit.timeit(__npdot, number=number)
+    assert_(time_dot < time_npdot)
+
+
+def test_cross():
+    # Test that cross is faster and equal to numpy.cross.
+    out = np.zeros(3)
+    for _ in range(10):
+        vec1 = np.random.random(3)
+        vec2 = np.random.random(3)
+        cross(out, vec1, vec2)
+        assert_array_equal(out, np.cross(vec1, vec2))
+
+    vec1 = np.random.random(3)
+    vec2 = np.random.random(3)
+
+    def __cross():
+        cross(out, vec1, vec2)
+
+    def __npcross():
+        np.cross(vec1, vec2)
+
+    number=10000
+    time_cross = timeit.timeit(__cross, number=number)
+    time_npcross = timeit.timeit(__npcross, number=number)
+    assert_(time_cross < time_npcross)