'''
pymaid_fiber_data.py
extract the muscle/tendon property from CATMAID data
suitable for muscles with associated tendons

input:
Name of the muscle in CATMAID project to extract data from
Name of the tendon in CATMAID project associated with the muscle
Name of the file to store the extracted data

output:
A csv file to store the information extracted in following format:

column_name:
['skid','fiber_name','fiber_len_integrated','fiber_len_convexHull','fiber_len_mainChain','fiber_tendon_angle','fiber_Zaxis_angle']
skid: skeleton id as defined in CATMAID
fiber_name: muscle fiber name inside a muscle cluster (one muscle cluster has multiple muscle fibers)
fiber_len_integrated: fiber length obtained with integration method by summing up length of each small segment
fiber_len_convexHull: fiber length obtained with furthest distance in between nodes in a convex hull wraping up the muscle (most general one to choose)
fiber_len_mainChain: fiber length obtained as distance between start and end node
fiber_tendon_angle: pennation angle in between individual muscle fibers with the corresponding tendon designated [for ACSA]
fiber_Zaxis_angle: pennation angle in between the muscle fiber and Z axis (orthorgnol to cross-section images) [for PCSA]

# requirement

pip install python-catmaid

# alert

# usage
get tibia MTU data

python pymaid_fiber_data.py -m "tibia extensor muscle fiber" -t "tibia extensor tendon" -o "tibia extensor.csv"
python pymaid_fiber_data.py -m "tibia flexor muscle fiber" -t "tibia flexor tendon" -o "tibia flexor.csv"

python pymaid_fiber_data_notendon.py -m "tibia extensor tendon"  -o "tibia extensor tendon.csv"
python pymaid_fiber_data_notendon.py -m "tibia flexor tendon"  -o "tibia flexor tendon.csv"

get femur MTU data

python pymaid_fiber_data.py -m "trochanter extensor muscle" -t "trochanter extensor tendon" -o "trochanter extensor muscle.csv"
python pymaid_fiber_data.py -m "sternotrochanter extensor muscle" -t "trochanter extensor tendon"  -o "sternotrochanter extensor muscle.csv"
python pymaid_fiber_data.py -m "tergotrochanteral muscle fiber" -t "trochanter extensor tendon"  -o "tergotrochanteral extensor muscle.csv"

python pymaid_fiber_data_notendon.py -m "trochanter flexor muscle"  -o "trochanter flexor muscle.csv"
python pymaid_fiber_data_notendon.py -m "accessory trochanter flexor muscle"  -o "accessory trochanter flexor muscle.csv"

python pymaid_fiber_data_notendon.py -m "trochanter extensor tendon"  -o "trochanter extensor tendon.csv"
python pymaid_fiber_data_notendon.py -m "trochanter flexor tendon"  -o "trochanter flexor tendon.csv"

get coxa fiber data

python pymaid_fiber_data_notendon.py -m "Pleural remotor and abductor (Miller 29)"  -o "Pleural remotor and abductor.csv"
python pymaid_fiber_data_notendon.py -m "pleural promotor (Miller 30)"  -o "pleural promotor.csv"
python pymaid_fiber_data_notendon.py -m "sternal posterior rotator (Miller 32)"  -o "sternal posterior rotator.csv"
python pymaid_fiber_data_notendon.py -m "sternal anterior rotator (Miller 31)"  -o "sternal anterior rotator.csv"
python pymaid_fiber_data_notendon.py -m "sternal adductor (Miller 33)"  -o "sternal abductor.csv"


'''
import pymaid
import numpy as np
import pandas as pd
from scipy import spatial
import argparse

## project and API keys for the project
##public project
pymaid_instance = pymaid.CatmaidInstance('https://radagast.hms.harvard.edu/catmaidvnc',
                                'd3e40f9bc75b4d70b558b1c75643e5b66e3342c0')
pymaid_instance.project_id = 61
# pymaid_instance = ##private project
# pymaid_instance = pymaid.CatmaidInstance(PROJECT_KEY)
# pymaid_instance.project_id = 30
## trochanter flexor tendon
## trochanter flexor muscle
## trochanter flexor.csv
parser = argparse.ArgumentParser(description='extract catmaid data into csv file')
parser.add_argument('-m', '--muscle')
parser.add_argument('-t', '--tendon')
parser.add_argument('-o', '--output')
args = parser.parse_args()

if args.muscle:
    muscle_name = args.muscle
else:
    muscle_name = 'tibia extensor muscle fiber'
print('Muscle: ', args.muscle)

if args.tendon:
    tendon_name = args.tendon
else:
    tendon_name = 'tibia extensor tendon'
print('Tendon: ', args.muscle)

if args.output:
    fout = args.output
else:
    fout = 'tibia extensor.csv'

print('Exporting muscle tendon data to ', fout)



def getZAngle(p1_start, p1_end):
    '''calculate fiber angle with Zaxis-plane, return degree'''
    ''' p_start: start point, p_end: end point'''
    # length_c
    length_c = get_p1p2_distance(p1_start,p1_end)
    # length_b
    p1_start_xy = p1_start
    p1_end_xy = p1_end
    p1_start_xy.z = 0
    p1_end_xy.z = 0
    length_b = get_p1p2_distance(p1_start_xy,p1_end_xy)
    cosAngle = length_b/length_c
    radian = np.arccos(cosAngle)
    return 90 - 180 / np.pi * radian

def getAngle(p1_start, p1_end, p2_start, p2_end):
    '''calculate pannetion vectors angle, return degree'''
    ''' p_start: start point, p_end: end point'''
    # vector 1
    x1 = list(p1_start.x)[0] - list(p1_end.x)[0]
    y1 = list(p1_start.y)[0] - list(p1_end.y)[0]
    z1 = list(p1_start.z)[0] - list(p1_end.z)[0]
    # vector 2
    x2 = list(p2_start.x)[0] - list(p2_end.x)[0]
    y2 = list(p2_start.y)[0] - list(p2_end.y)[0]
    z2 = list(p2_start.z)[0] - list(p2_end.z)[0]
    vectorDot = x1*x2 + y1*y2 + z1*z2
    vectorMold1 = np.sqrt(x1*x1 + y1*y1 + z1*z1)
    vectorMold2 = np.sqrt(x2*x2 + y2*y2 + z2*z2)
    cosAngle = vectorDot / (vectorMold1 * vectorMold2)
    radian = np.arccos(cosAngle)
    return 180 / np.pi * radian

def get_p1p2_distance(p1,p2):
    '''get distance with p1 and p2'''
    x1 = p1.iloc[0]['x']
    y1 = p1.iloc[0]['y']
    z1 = p1.iloc[0]['z']
    x2 = p2.iloc[0]['x']
    y2 = p2.iloc[0]['y']
    z2 = p2.iloc[0]['z']
    distance = np.sqrt((x1-x2)**2 + (y1-y2)**2 + (z1-z2)**2)
    return(distance)

def get_distance_convexHull(tn_table):
    '''get two points from all points tn_table, ConvexHull method'''
    M = tn_table[['x','y','z','skeleton_id']]
    pts = tn_table[['x','y','z']]
    # two points which are fruthest apart will occur as vertices of the convex hull
    idx = spatial.ConvexHull(pts).vertices
    candidates = pts.T[idx].T
    # get distances between each pair of candidate points
    dist_mat = spatial.distance_matrix(candidates,candidates)
    # get indices of candidates that are furthest apart
    i,j = np.unravel_index(dist_mat.argmax(),dist_mat.shape)
    p1_idx = candidates.index[i]
    p2_idx = candidates.index[j]
    p1 = M[p1_idx:p1_idx+1]
    p2 = M[p2_idx:p2_idx+1]
    M1 = M.T[[p1_idx,p2_idx]].T[['x','y','z']]
    dist_mat1 = spatial.distance_matrix(M1,M1)
    distance = dist_mat1[0][1]
    return(p1,p2,distance)

def get_distance_integrated(fiber_table):
    '''get integrated length for fiber '''
    # not continued   skid=609306 segment node_id=12592618 has no child,and node_id 12592524 has two children(12592619,12592525)
    fiber_length = 0.0
    # p2='node_id', p1='parent_node_id'
    for i in range(len(fiber_table)):
        node_id_p1 = fiber_table.iloc[i]['parent_node_id']
        if node_id_p1==-1:
            continue
        node_id_p2 = fiber_table.iloc[i]['node_id']
        p1 = fiber_table[fiber_table['node_id']==node_id_p1]
        p2 = fiber_table[fiber_table['node_id']==node_id_p2]
        segment_legth = get_p1p2_distance(p1,p2)
        fiber_length += segment_legth
    return(fiber_length)

def get_distance_mainChain(fiber_table):
    '''get length of start_point and end_point in the main link_chain'''
    # 'node_id', 'parent_node_id'  segment_count = len - 1, start with pid=-1
    node_id_p1 = -1
    p0 = fiber_table[fiber_table['parent_node_id']==node_id_p1]
    # find tail
    p1 = p0
    for i in range(len(fiber_table)-1):
        node_id_p1 = p1.iloc[0]['node_id']
        p2 = fiber_table[fiber_table['parent_node_id']==node_id_p1]
        if len(p2)!=0:
            p1 = p2
        else:
            exit
    # print("found main_chain's tail:{}\n".format(node_id_p1))
    fiber_length = get_p1p2_distance(p0,p1)
    return(fiber_length)


tendon_skids = pymaid.get_skids_by_annotation(tendon_name)
table = pymaid.get_node_table(tendon_skids)
tendon_length = get_distance_integrated(table)
print('distance_integrated:',tendon_length)
tendon_length = get_distance_mainChain(table)
print('distance_mainChain:',tendon_length)
tendon_start,tendon_end,distance = get_distance_convexHull(table)
print('distance_convexHull:{}\n'.format(distance))
print('tendon_start_point is :\n',tendon_start)
print('tendon_start_point is :\n',tendon_end)

skids = pymaid.get_skids_by_annotation(muscle_name)
fibers = []
for skid in skids:
    try:
        fiber_table = pymaid.get_node_table(skid)
        fiber_len_integrated = get_distance_integrated(fiber_table)
        fiber_len_mainChain = get_distance_mainChain(fiber_table)
        p_start,p_end,fiber_len_convexHull = get_distance_convexHull(fiber_table)
        fiber_tendon_angle = getAngle(p_start, p_end, tendon_start, tendon_end)
        fiber_Zaxis_angle = getZAngle(p_start, p_end)
        fiber_name = pymaid.get_names(skid).get(str(skid))
        fibers.append([skid,fiber_name,fiber_len_integrated,fiber_len_convexHull,
            fiber_len_mainChain,fiber_tendon_angle,fiber_Zaxis_angle])
    except Exception as exception:
        print(f'Unable to create tree node table for skid {skid} due to {type(exception).__name__}, skipping...')
fibers_df = pd.DataFrame(fibers,columns=['skid','fiber_name','fiber_len_integrated',
                               'fiber_len_convexHull','fiber_len_mainChain','fiber_tendon_angle','fiber_Zaxis_angle'])
fibers_df.to_csv(fout)
print("fibers' data has been write to file ", args.output)

exit()



