''' Metrics for tracks, where tracks are arrays of points '''
from __future__ import division, print_function, absolute_import

from ..utils.six.moves import xrange

import numpy as np
from scipy.interpolate import splprep, splev


def winding(xyz):
    '''Total turning angle projected.

    Project space curve to best fitting plane. Calculate the cumulative signed
    angle between each line segment and the previous one.

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

    Returns
    ---------
    a : scalar
        Total turning angle in degrees.

    '''

    U, s, V = np.linalg.svd(xyz-np.mean(xyz, axis=0), 0)
    proj = np.dot(U[:, 0:2], np.diag(s[0:2]))

    turn = 0
    for j in range(len(xyz)-1):
        v0 = proj[j]
        v1 = proj[j+1]
        v = np.dot(v0, v1)/(np.linalg.norm(v0)*np.linalg.norm(v1))
        v = np.clip(v, -1, 1)
        tmp = np.arccos(v)
        turn += tmp

    return np.rad2deg(turn)


def length(xyz, along=False):
    ''' Euclidean length of track line

    This will give length in mm if tracks are expressed in world coordinates.

    Parameters
    ------------
    xyz : array-like shape (N,3)
       array representing x,y,z of N points in a track
    along : bool, optional
       If True, return array giving cumulative length along track,
       otherwise (default) return scalar giving total length.

    Returns
    ---------
    L : scalar or array shape (N-1,)
       scalar in case of `along` == False, giving total length, array if
       `along` == True, giving cumulative lengths.

    Examples
    ----------
    >>> from dipy.tracking.metrics import length
    >>> xyz = np.array([[1,1,1],[2,3,4],[0,0,0]])
    >>> expected_lens = np.sqrt([1+2**2+3**2, 2**2+3**2+4**2])
    >>> length(xyz) == expected_lens.sum()
    True
    >>> len_along = length(xyz, along=True)
    >>> np.allclose(len_along, expected_lens.cumsum())
    True
    >>> length([])
    0
    >>> length([[1, 2, 3]])
    0
    >>> length([], along=True)
    array([0])
    '''
    xyz = np.asarray(xyz)
    if xyz.shape[0] < 2:
        if along:
            return np.array([0])
        return 0
    dists = np.sqrt((np.diff(xyz, axis=0)**2).sum(axis=1))
    if along:
        return np.cumsum(dists)
    return np.sum(dists)


def bytes(xyz):
    '''Size of track in bytes.

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

    Returns
    ---------
    b : int
        Number of bytes.

    '''
    return xyz.nbytes


def midpoint(xyz):
    ''' Midpoint of track

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

    Returns
    ---------
    mp : array shape (3,)
       Middle point of line, such that, if L is the line length then
       `np` is the point such that the length xyz[0] to `mp` and from
       `mp` to xyz[-1] is L/2.  If the middle point is not a point in
       `xyz`, then we take the interpolation between the two nearest
       `xyz` points.  If `xyz` is empty, return a ValueError

    Examples
    --------
    >>> from dipy.tracking.metrics import midpoint
    >>> midpoint([])
    Traceback (most recent call last):
       ...
    ValueError: xyz array cannot be empty
    >>> midpoint([[1, 2, 3]])
    array([1, 2, 3])
    >>> xyz = np.array([[1,1,1],[2,3,4]])
    >>> midpoint(xyz)
    array([ 1.5,  2. ,  2.5])
    >>> xyz = np.array([[0,0,0],[1,1,1],[2,2,2]])
    >>> midpoint(xyz)
    array([ 1.,  1.,  1.])
    >>> xyz = np.array([[0,0,0],[1,0,0],[3,0,0]])
    >>> midpoint(xyz)
    array([ 1.5,  0. ,  0. ])
    >>> xyz = np.array([[0,9,7],[1,9,7],[3,9,7]])
    >>> midpoint(xyz)
    array([ 1.5,  9. ,  7. ])
    '''
    xyz = np.asarray(xyz)
    n_pts = xyz.shape[0]
    if n_pts == 0:
        raise ValueError('xyz array cannot be empty')
    if n_pts == 1:
        return xyz.copy().squeeze()
    cumlen = np.zeros(n_pts)
    cumlen[1:] = length(xyz, along=True)
    midlen = cumlen[-1]/2.0
    ind = np.where((cumlen-midlen) > 0)[0][0]
    len0 = cumlen[ind-1]
    len1 = cumlen[ind]
    Ds = midlen-len0
    Lambda = Ds/(len1-len0)
    return Lambda*xyz[ind]+(1-Lambda)*xyz[ind-1]


def center_of_mass(xyz):
    ''' Center of mass of streamline

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

    Returns
    ---------
    com : array shape (3,)
       center of mass of streamline

    Examples
    ----------
    >>> from dipy.tracking.metrics import center_of_mass
    >>> center_of_mass([])
    Traceback (most recent call last):
       ...
    ValueError: xyz array cannot be empty
    >>> center_of_mass([[1,1,1]])
    array([ 1.,  1.,  1.])
    >>> xyz = np.array([[0,0,0],[1,1,1],[2,2,2]])
    >>> center_of_mass(xyz)
    array([ 1.,  1.,  1.])
    '''
    xyz = np.asarray(xyz)
    if xyz.size == 0:
        raise ValueError('xyz array cannot be empty')
    return np.mean(xyz, axis=0)


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)


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

    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)

    Parameters
    ----------
    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