def frenet_serret(xyz):
    ''' Frenet-Serret Space Curve Invariants

    Taken from "Diffusion Imaging in Python" `DiPy package<http://nipy.org/dipy/>`_

    Calculates the 3 vector and 2 scalar invariants of a space curve defined by vectors r = (x,y,z).  If z is omitted (i.e. the array xyz has shape (N,2), then the curve is only 2D (planar), but the equations are still valid.

    Similar to http://www.mathworks.com/matlabcentral/fileexchange/11169

    In the following equations the prime ($'$) indicates differentiation with respect to the parameter $s$ of a parametrised curve $\mathbf{r}(s)$.

    - $\mathbf{T}=\mathbf{r'}/|\mathbf{r'}|\qquad$ (Tangent vector)}

    - $\mathbf{N}=\mathbf{T'}/|\mathbf{T'}|\qquad$ (Normal vector)

    - $\mathbf{B}=\mathbf{T}\times\mathbf{N}\qquad$ (Binormal vector)

    - $\kappa=|\mathbf{T'}|\qquad$ (Curvature)

    - $\mathrm{\tau}=-\mathbf{B'}\cdot\mathbf{N}$ (Torsion)

    **Arguments:**

            * xyz : array-like shape (N,3)
                            array representing x,y,z of N points in a track

    **Returns:**

                    * T : array shape (N,3)
                            array representing the tangent of the curve xyz
                    * N : array shape (N,3)
                            array representing the normal of the curve xyz
            * B : array shape (N,3)
                    array representing the binormal of the curve xyz
            * k : array shape (N,1)
                    array representing the curvature of the curve xyz
            * t : array shape (N,1)
                    array representing the torsion of the curve xyz

    Examples:

    Create a helix and calculate its tangent, normal, binormal, curvature and torsion

    >>> from dipy.tracking import metrics as tm
    >>> import numpy as np
    >>> theta = 2*np.pi*np.linspace(0,2,100)
    >>> x=np.cos(theta)
    >>> y=np.sin(theta)
    >>> z=theta/(2*np.pi)
    >>> xyz=np.vstack((x,y,z)).T
    >>> T,N,B,k,t=tm.frenet_serret(xyz)

    '''

    def magn(xyz, n=1):
        ''' magnitude of vector


        '''
        mag = np.sum(xyz**2, axis=1)**0.5
        imag = np.where(mag == 0)
        mag[imag] = np.finfo(float).eps

        if n > 1:
            return np.tile(mag, (n, 1)).T

        return mag.reshape(len(mag), 1)

    xyz = np.asarray(xyz)
    n_pts = xyz.shape[0]

    if n_pts == 0:
        raise ValueError('xyz array cannot be empty')

    dxyz = np.gradient(xyz)[0]
    ddxyz = np.gradient(dxyz)[0]

    # Tangent
    T = np.divide(dxyz, magn(dxyz, 3))

    # Derivative of Tangent
    dT = np.gradient(T)[0]

    # Normal
    N = np.divide(dT, magn(dT, 3))

    # Binormal
    B = np.cross(T, N)

    # Curvature
    k = magn(np.cross(dxyz, ddxyz), 1) / (magn(dxyz, 1)**3)

    # Torsion
    #(In matlab was t=dot(-B,N,2))
    t = np.sum(-B * N, axis=1)
    # return T,N,B,k,t,dxyz,ddxyz,dT

    return T, N, B, k, t
Frenet-Serret Space Curve Invariants

    Taken from "Diffusion Imaging in Python" `DiPy package<http://nipy.org/dipy/>`_

    Calculates the 3 vector and 2 scalar invariants of a space curve defined by vectors r = (x,y,z).  If z is omitted (i.e. the array xyz has shape (N,2), then the curve is only 2D (planar), but the equations are still valid.

    Similar to http://www.mathworks.com/matlabcentral/fileexchange/11169

    In the following equations the prime ($'$) indicates differentiation with respect to the parameter $s$ of a parametrised curve $\mathbf{r}(s)$.

    - $\mathbf{T}=\mathbf{r'}/|\mathbf{r'}|\qquad$ (Tangent vector)}

    - $\mathbf{N}=\mathbf{T'}/|\mathbf{T'}|\qquad$ (Normal vector)

    - $\mathbf{B}=\mathbf{T}\times\mathbf{N}\qquad$ (Binormal vector)

    - $\kappa=|\mathbf{T'}|\qquad$ (Curvature)

    - $\mathrm{\tau}=-\mathbf{B'}\cdot\mathbf{N}$ (Torsion)

    **Arguments:**

            * xyz : array-like shape (N,3)
                            array representing x,y,z of N points in a track

    **Returns:**

                    * T : array shape (N,3)
                            array representing the tangent of the curve xyz
                    * N : array shape (N,3)
                            array representing the normal of the curve xyz
            * B : array shape (N,3)
                    array representing the binormal of the curve xyz
            * k : array shape (N,1)
                    array representing the curvature of the curve xyz
            * t : array shape (N,1)
                    array representing the torsion of the curve xyz

    Examples:

    Create a helix and calculate its tangent, normal, binormal, curvature and torsion

    >>> from dipy.tracking import metrics as tm
    >>> import numpy as np
    >>> theta = 2*np.pi*np.linspace(0,2,100)
    >>> x=np.cos(theta)
    >>> y=np.sin(theta)
    >>> z=theta/(2*np.pi)
    >>> xyz=np.vstack((x,y,z)).T
    >>> T,N,B,k,t=tm.frenet_serret(xyz)