In [1]:
import numpy as np
import scipy.linalg
import matplotlib
import matplotlib.pyplot as plt
from scipy.io import loadmat
import pandas as pd
import seaborn as sns
%matplotlib inline
matplotlib.rcParams.update({'font.size': 14})
In [2]:
skeleton_data=loadmat('data/skeletons_OPL_final.mat')
skeleton_ids=skeleton_data['kn_allSkeletons_e2006_IDs'].flatten()
skeletons=skeleton_data['kn_allSkeletons_e2006'].flatten()
M=np.loadtxt('data/coordinate_rotation')
BC_ids=np.loadtxt('data/BC_IDs_new').astype(int)
In [3]:
nodes_complete=[]
edges_complete=[]
for i in range(skeleton_ids[skeleton_ids<2000].shape[0]):
nodes=skeletons[i].item()[list(skeletons[i].dtype.names).index('nodes')][:,:3]
edges=skeletons[i].item()[list(skeletons[i].dtype.names).index('edges')]
nodes=nodes[np.unique(edges)-1,:]
edge_numbers=np.unique(edges)
edge_mapping=np.vstack((edge_numbers,np.arange(1,edge_numbers.shape[0]+1))).T
edges_new=np.zeros(edges.shape)
for k in range(edges.shape[0]):
edges_new[k,0]=edge_mapping[edge_mapping[:,0]==edges[k,0],1]
edges_new[k,1]=edge_mapping[edge_mapping[:,0]==edges[k,1],1]
nodes_complete.append(nodes)
edges_complete.append(edges_new)
In [4]:
contact_data=loadmat('data/contacts_final.mat')['kn_allContactData_Interfaces']
cone_centers=np.loadtxt('data/cone_means_rot.txt')
cone_boundaries=np.loadtxt('data/cone_boundaries_rot.txt')
cone_boundaries[:,1:]*=16.5
Columns in contact_data_cones:
0: Cell ID
1: Cone ID
2: Type ID
3-5: Contact position (absolute)
6-8: Contact position relative to cone center
9: Excentricity (radial distance to cone center)
10: Distance to cone top
11: Distance to cone bottom
12: Position relative to cone boundaries
13: Contact area
14: Branching distance
15: Ending distance
16: Angle
In [5]:
contact_data_red=contact_data[(contact_data[:,1]>2000)&(contact_data[:,1]<3000)&np.in1d(contact_data[:,0],BC_ids[:,0])]
contact_data_red=contact_data_red[contact_data_red[:,4]>64,:]
contact_data_red=contact_data_red[contact_data_red[:,4]<192,:]
contact_data_red=contact_data_red[contact_data_red[:,5]>64,:]
contact_data_red=contact_data_red[contact_data_red[:,5]<192,:]
contact_data_red=contact_data_red[contact_data_red[:,6]>64,:]
contact_data_red=contact_data_red[contact_data_red[:,6]<192,:]
In [6]:
#set up and populate area
contact_data_cones=np.hstack((contact_data_red[:,:2],np.zeros((contact_data_red.shape[0],1)),np.arange(contact_data_red.shape[0]).reshape(-1,1),\
(contact_data_red[:,4]+contact_data_red[:,11]+64).reshape(-1,1),\
(contact_data_red[:,5]+contact_data_red[:,12]+64).reshape(-1,1),(contact_data_red[:,6]+contact_data_red[:,13]+64).reshape(-1,1),\
np.zeros((contact_data_red.shape[0],7)),contact_data_red[:,14].reshape(-1,1),np.zeros((contact_data_red.shape[0],3))))
In [7]:
#populate type ids (new classification)
for i in range(contact_data_cones.shape[0]):
contact_data_cones[i,2]=BC_ids[np.where(BC_ids[:,0]==contact_data_cones[i,0]),:].ravel()[4]
In [8]:
#populate branching and ending distances - has to be done before coordinate rotation!
for BC in np.unique(contact_data_cones[:,0]):
nodes=[]
edges=[]
for cell in np.where(skeleton_ids==BC)[0]:
nodes.append(nodes_complete[cell])
edges.append(edges_complete[cell])
num_skel=len(nodes)
for i in np.where(contact_data_cones[:,0]==BC)[0]:
contact=contact_data_cones[i,4:7]*[16.5,16.5,25]
branching_distances=[]
ending_distances=[]
for skel in range(num_skel):
contact_node_dist=np.linalg.norm(nodes[skel][:,:3]-contact,axis=1)
if np.min(contact_node_dist)>1000:
continue
contact_node=np.where(contact_node_dist==np.min(contact_node_dist))[0][0]+1
distance_to_branching=[]
distance_to_ending=[]
dist=0
nodes_captured=[contact_node]
node=contact_node
while True:
next_nodes=np.unique(edges[skel][np.where(edges[skel][:,:]==node)[0],:]).astype(int)
if next_nodes.shape[0]>=4:
distance_to_branching.append(dist)
break
elif next_nodes.shape[0]==2:
distance_to_ending.append(dist)
if not node==contact_node:
break
next_nodes=np.setdiff1d(next_nodes,nodes_captured)
if next_nodes.shape[0]==0:
break
next_node=next_nodes[0]
dist+=np.linalg.norm(nodes[skel][node-1,:3]-nodes[skel][next_node-1,:3])
nodes_captured.append(next_node)
node=next_node
dist=0
node=contact_node
while True:
next_nodes=np.unique(edges[skel][np.where(edges[skel][:,:]==node)[0],:]).astype(int)
if next_nodes.shape[0]>=4:
distance_to_branching.append(dist)
break
elif next_nodes.shape[0]==2:
distance_to_ending.append(dist)
break
next_nodes=np.setdiff1d(next_nodes,nodes_captured)
if next_nodes.shape[0]==0:
break
next_node=next_nodes[0]
dist+=np.linalg.norm(nodes[skel][node-1,:3]-nodes[skel][next_node-1,:3])
nodes_captured.append(next_node)
node=next_node
if not len(distance_to_branching)==0:
branching_distances.append(min(distance_to_branching))
if not len(distance_to_ending)==0:
ending_distances.append(min(distance_to_ending))
if not len(branching_distances)==0:
contact_data_cones[i,15]=np.mean(branching_distances)
else:
contact_data_cones[i,15]=100000
if not len(ending_distances)==0:
contact_data_cones[i,16]=np.mean(ending_distances)
else:
contact_data_cones[i,16]=100000
In [9]:
#populate angle
contact_node_edges=[]
for BC in np.unique(contact_data_cones[:,0]):
nodes=[]
edges=[]
for cell in np.where(skeleton_ids==BC)[0]:
nodes.append(nodes_complete[cell])
edges.append(edges_complete[cell])
num_skel=len(nodes)
for i in np.where(contact_data_cones[:,0]==BC)[0]:
contact=contact_data_cones[i,4:7]*[16.5,16.5,25]
for skel in range(num_skel):
contact_node_dist=np.linalg.norm(nodes[skel][:,:3]-contact,axis=1)
if np.min(contact_node_dist)>1000:#300:
continue
contact_node=np.where(contact_node_dist==np.min(contact_node_dist))[0][0]+1
dist=0
nodes_captured=[contact_node]
node=contact_node
test_nodes=np.unique(edges[skel][(edges[skel][:,0]==node)|(edges[skel][:,1]==node),:])
test_nodes=test_nodes[test_nodes!=node].astype(int)
test_dist=np.linalg.norm(nodes[skel][test_nodes-1,:3]-contact,axis=1)
node2=test_nodes[test_dist==np.min(test_dist)]
test_nodes=[contact_node-1,node2[0]-1]
rel_nodes=nodes[skel][test_nodes,:3]/[16.5,16.5,25]
rel_nodes=np.dot(M,rel_nodes.T).T*[16.5,16.5,25]
rel_vec=rel_nodes[0,:]-rel_nodes[1,:]
contact_data_cones[i,17]=np.arctan(np.abs(rel_vec[0])/np.sqrt(rel_vec[1]**2+rel_vec[2]**2))
In [10]:
#rotate coordinates
contact_data_cones[:,4:7]=np.dot(M,contact_data_cones[:,4:7].T).T
for i in range(contact_data_cones.shape[0]):
contact_data_cones[i,4:7]=contact_data_cones[i,4:7]*[16.5,16.5,25]
In [11]:
#populate position relative to cone center
for i in range(contact_data_cones.shape[0]):
contact_data_cones[i,7:10]=contact_data_cones[i,4:7]-cone_centers[cone_centers[:,0]==contact_data_cones[i,1],1:4].ravel()*[16.5,16.5,25]
In [12]:
#populate excentricity
contact_data_cones[:,10]=np.sqrt(contact_data_cones[:,8]**2+contact_data_cones[:,9]**2)
In [13]:
#populate distance to top and bottom and height
for i in range(contact_data_cones.shape[0]):
boundaries=cone_boundaries[cone_boundaries[:,0]==contact_data_cones[i,1],1:].flatten()
contact_data_cones[i,11]=contact_data_cones[i,4]-boundaries[0]
contact_data_cones[i,12]=contact_data_cones[i,4]-boundaries[1]
contact_data_cones[i,13]=contact_data_cones[i,12]/(boundaries[0]-boundaries[1])
In [14]:
all_contacts=pd.DataFrame(contact_data_cones,columns=['cell','cone','type','contact_id','x_pos','y_pos','z_pos',\
'x_center_distance','y_center_distance','z_center_distance','excentricity',\
'distance_to_upper_boundary','distance_to_lower_boundary','height','area','branch_dist','end_dist','angle'])
In [ ]:
# all_contacts.to_pickle('../data/cone_contact_data')
In [ ]: