Rod contact analysis

This notebook contains the code to collect the data used for the contact classification of rod to bipolar cell contacts



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'][0]
M=np.loadtxt('data/coordinate_rotation')
BC_ids=np.loadtxt('data/BC_IDs_new')



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']



In [5]:

    
rod_centers=np.loadtxt('data/rod_means_rot.txt')



In [6]:

    
rod_boundaries=np.loadtxt('data/rod_boundaries_rot.txt')
rod_boundaries[:,1:]*=16.5

Columns in contact_data_rods:

0: Cell ID
1: Rod ID
2: Type ID
3-5: Contact position (absolute)
6-8: Contact position relative to rod center
9: Excentricity (radial distance to rod center)
10: Distance to rod top
11: Distance to rod bottom
12: Position relative to rod boundaries
13: Contact area
14: Branching distance
15: Ending distance
16: Angle



In [7]:

    
contact_data_red=contact_data[(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 [8]:

    
#set up and populate area
contact_data_rods=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 [9]:

    
#populate type ids (new classification)
for i in range(contact_data_rods.shape[0]):
    contact_data_rods[i,2]=BC_ids[np.where(BC_ids[:,0]==contact_data_rods[i,0]),:].ravel()[4]



In [10]:

    
#populate branching and ending distances - has to be done before coordinate rotation!
for BC in np.unique(contact_data_rods[:,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_rods[:,0]==BC)[0]:
        contact=contact_data_rods[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_rods[i,15]=np.mean(branching_distances)
        else:
            contact_data_rods[i,15]=100000
        if not len(ending_distances)==0:
            contact_data_rods[i,16]=np.mean(ending_distances)
        else:
            contact_data_rods[i,16]=100000



In [11]:

    
#populate angle
contact_node_edges=[]
for BC in np.unique(contact_data_rods[:,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_rods[:,0]==BC)[0]:
        contact=contact_data_rods[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_rods[i,17]=np.arctan(np.abs(rel_vec[0])/np.sqrt(rel_vec[1]**2+rel_vec[2]**2))



In [12]:

    
#rotate coordinates
contact_data_rods[:,4:7]=np.dot(M,contact_data_rods[:,4:7].T).T
for i in range(contact_data_rods.shape[0]):
    contact_data_rods[i,4:7]=contact_data_rods[i,4:7]*[16.5,16.5,25]



In [13]:

    
#populate position relative to rod center
for i in range(contact_data_rods.shape[0]):
    contact_data_rods[i,7:10]=contact_data_rods[i,4:7]-rod_centers[rod_centers[:,0]==contact_data_rods[i,1],1:4].ravel()*[16.5,16.5,25]



In [14]:

    
#populate excentricity
contact_data_rods[:,10]=np.sqrt(contact_data_rods[:,8]**2+contact_data_rods[:,9]**2)



In [15]:

    
#populate distance to top and bottom and height
for i in range(contact_data_rods.shape[0]):
    boundaries=rod_boundaries[rod_boundaries[:,0]==contact_data_rods[i,1],1:].flatten()
    contact_data_rods[i,11]=contact_data_rods[i,4]-boundaries[0]
    contact_data_rods[i,12]=contact_data_rods[i,4]-boundaries[1]
    contact_data_rods[i,13]=contact_data_rods[i,12]/(boundaries[0]-boundaries[1])



In [16]:

    
all_contacts=pd.DataFrame(contact_data_rods,columns=['cell','rod','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 [17]:

    
# all_contacts.to_pickle('data/rod_contact_data')

Identify bipolar cells within the area of the reconstructed rods



In [18]:

    
BC_means=[]
for cell in np.unique(skeleton_ids[skeleton_ids<2000]):
    mean=[]
    for i in np.where(skeleton_ids==cell)[0]:
        nodes=skeletons[i].item()[list(skeletons[i].dtype.names).index('nodes')][:,:3]/[16.5,16.5,25]
        nodes=nodes[np.dot(M,nodes.T).T[:,0]<4000]
        if nodes.shape[0]>1:
            mean.append(np.mean(nodes[:,1:],axis=0))
    mean=np.array(mean).reshape(-1,2)
    if mean.shape[0]>0:
        BC_means.append([cell,np.mean(mean[:,0]),np.mean(mean[:,1])])
BC_means=np.array(BC_means)



In [19]:

    
rod_area_max=np.zeros(2)
rod_area_min=np.ones(2)*10000
for i in range(skeleton_ids.shape[0]):
    if skeleton_ids[i]>3000:
        nodes=skeletons[i].item()[list(skeletons[i].dtype.names).index('nodes')][:,:3]/[16.5,16.5,25]
        
        if nodes.shape[0]>0:
            rod_area_max=np.maximum(rod_area_max,np.max(nodes[:,1:],axis=0))
            rod_area_min=np.minimum(rod_area_min,np.min(nodes[:,1:],axis=0))



In [20]:

    
BC_in_rod_area=BC_means[np.all(BC_means[:,1:]>rod_area_min,axis=1)&np.all(BC_means[:,1:]<rod_area_max,axis=1),0]



In [21]:

    
# np.savetxt('data/BC_in_rod_area',BC_in_rod_area)



In [ ]: