cone_contact_analysis.ipynb

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