'''
Utilities for analyzing movement data recorded with marker-based video
systems.
'''

'''
Author: Thomas Haslwanter
Version: 1.1
Date: April-2016
'''

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import pandas as pd 
from numpy import r_, sum

# The following construct is required since I want to run the module as a script
# inside the thLib-directory
import os
import sys
PACKAGE_PARENT = '..'
SCRIPT_DIR = os.path.realpath(os.path.dirname(__file__))
sys.path.append(os.path.normpath(os.path.join(SCRIPT_DIR, PACKAGE_PARENT)))

from skinematics import quat, vector

def analyze3Dmarkers(MarkerPos, ReferencePos):
    '''
    Take recorded positions from 3 markers, and calculate center-of-mass (COM) and orientation
    Can be used e.g. for the analysis of Optotrac data.

    Parameters
    ----------
    MarkerPos : ndarray, shape (N,9)
        x/y/z coordinates of 3 markers

    ReferencePos : ndarray, shape (1,9)
        x/y/z coordinates of markers in the reference position

    Returns
    -------
    Position : ndarray, shape (N,3)
        x/y/z coordinates of COM, relative to the reference position
    Orientation : ndarray, shape (N,3)
        Orientation relative to reference orientation, expressed as quaternion

    Example
    -------
    >>> (PosOut, OrientOut) = analyze3Dmarkers(MarkerPos, ReferencePos)


    '''

    # Specify where the x-, y-, and z-position are located in the data
    cols = {'x' : r_[(0,3,6)]} 
    cols['y'] = cols['x'] + 1
    cols['z'] = cols['x'] + 2    

    # Calculate the position
    cog = np.vstack(( sum(MarkerPos[:,cols['x']], axis=1),
                      sum(MarkerPos[:,cols['y']], axis=1),
                      sum(MarkerPos[:,cols['z']], axis=1) )).T/3.

    cog_ref = np.vstack(( sum(ReferencePos[cols['x']]),
                          sum(ReferencePos[cols['y']]),
                          sum(ReferencePos[cols['z']]) )).T/3.                      

    position = cog - cog_ref    

    # Calculate the orientation    
    numPoints = len(MarkerPos)
    orientation = np.ones((numPoints,3))

    refOrientation = vector.plane_orientation(ReferencePos[:3], ReferencePos[3:6], ReferencePos[6:])

    for ii in range(numPoints):
        '''The three points define a triangle. The first rotation is such
        that the orientation of the reference-triangle, defined as the
        direction perpendicular to the triangle, is rotated along the shortest
        path to the current orientation.
        In other words, this is a rotation outside the plane spanned by the three
        marker points.'''

        curOrientation = vector.plane_orientation( MarkerPos[ii,:3], MarkerPos[ii,3:6], MarkerPos[ii,6:])
        alpha = vector.angle(refOrientation, curOrientation)        

        if alpha > 0:
            n1 = vector.normalize( np.cross(refOrientation, curOrientation) )
            q1 = n1 * np.sin(alpha/2)
        else:
            q1 = r_[0,0,0]

        # Now rotate the triangle into this orientation ...
        refPos_after_q1 = vector.rotate_vector(np.reshape(ReferencePos,(3,3)), q1)

        '''Find which further rotation in the plane spanned by the three marker points
	is necessary to bring the data into the measured orientation.'''

        Marker_0 = MarkerPos[ii,:3]
        Marker_1 = MarkerPos[ii,3:6]
        Vector10 = Marker_0 - Marker_1
        vector10_ref = refPos_after_q1[0]-refPos_after_q1[1]
        beta = vector.angle(Vector10, vector10_ref)

        q2 = curOrientation * np.sin(beta/2)

        if np.cross(vector10_ref,Vector10).dot(curOrientation)<=0:
            q2 = -q2
        orientation[ii,:] = quat.quatmult(q2, q1)

    return (position, orientation)

def findTrajectory(r0, Position, Orientation):
    '''
    Movement trajetory of a point on an object, from the position and
    orientation of a sensor, and the relative position of the point at t=0.

    Parameters
    ----------
    r0 : ndarray (3,)
        Position of point relative to center of markers, when the object is
        in the reference position.
    Position : ndarray, shape (N,3)
        x/y/z coordinates of COM, relative to the reference position
    Orientation : ndarray, shape (N,3)
        Orientation relative to reference orientation, expressed as quaternion

    Returns
    -------
    mov : ndarray, shape (N,3)
        x/y/z coordinates of the position on the object, relative to the
        reference position of the markers

    Notes
    ----- 

      .. math::

          \\vec C(t) = \\vec M(t) + \\vec r(t) = \\vec M(t) +
          {{\\bf{R}}_{mov}}(t) \\cdot \\vec r({t_0})

    Examples
    --------
    >>> t = np.arange(0,10,0.1)
    >>> translation = (np.c_[[1,1,0]]*t).T
    >>> M = np.empty((3,3))
    >>> M[0] = np.r_[0,0,0]
    >>> M[1]= np.r_[1,0,0]
    >>> M[2] = np.r_[1,1,0]
    >>> M -= np.mean(M, axis=0) 
    >>> q = np.vstack((np.zeros_like(t), np.zeros_like(t),quat.deg2quat(100*t))).T
    >>> M0 = vector.rotate_vector(M[0], q) + translation
    >>> M1 = vector.rotate_vector(M[1], q) + translation
    >>> M2 = vector.rotate_vector(M[2], q) + translation
    >>> data = np.hstack((M0,M1,M2))
    >>> (pos, ori) = signals.analyze3Dmarkers(data, data[0])
    >>> r0 = np.r_[1,2,3]
    >>> movement = movement_from_markers(r0, pos, ori)

    '''

    return Position + vector.rotate_vector(r0, Orientation)

if __name__ == '__main__':
    # Here a marker-test is required

    print('Done')
