cone_overview_projection.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cone projection images\n",
    "\n",
    "This notebook contains the code to draw all images showing projections of the cone volumes along the optical axis such as figures 1C&D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.io import loadmat\n",
    "from scipy.spatial import ConvexHull\n",
    "from PIL import Image\n",
    "from PIL import ImageDraw\n",
    "from PIL import ImageFont\n",
    "from shapely.geometry import Point\n",
    "from shapely.geometry import Polygon\n",
    "import pandas as pd\n",
    "%matplotlib inline\n",
    "matplotlib.rcParams.update({'font.size': 14})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "matplotlib.rcParams['font.family'] = 'sans-serif'\n",
    "matplotlib.rcParams['font.sans-serif'] = ['Arial']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "M=np.loadtxt('data/coordinate_rotation')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "soma_data=loadmat('../data/Soma_positions.mat')\n",
    "soma_pos=soma_data['kn_e2006_ALLSKELETONS_FINAL2012_allSomata']\n",
    "soma_internal_ids=soma_data['kn_e2006_ALLSKELETONS_FINAL2012_cellIDs']\n",
    "soma_line_ids=soma_data['kn_e2006_ALLSKELETONS_FINAL2012_cellIDs_pure_forSomata']\n",
    "soma_pos=np.dot(M,soma_pos[:,:3].T).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cut_off_cones=np.array([73,72,71,74,201,202,211,203,88,204,142,34,4,30,199])+2000\n",
    "cones_outside_cbc9=np.array([193,209,200,198,197,199,73,72,71,74,69,67,66,211,86,88,87,120,85,204,84,207,128,114,126,127,125,142,130])+2000\n",
    "Blue_cone_ids=np.array([2006,2007,2009,2024,2028,2029])\n",
    "Blue_cone_ids_alt=np.array([2007,2006,2010,2029,2061,2009,2031,2014,2024,2023,2100,2028,2017,2090])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "BC_ids=np.loadtxt('data/BC_IDs_new')\n",
    "\n",
    "cell_ids=[]\n",
    "for i in range(BC_ids.shape[0]):\n",
    "    if BC_ids[i,4] in range(58,72):\n",
    "        cell_ids.append([BC_ids[i,0],BC_ids[i,1],BC_ids[i,4]])\n",
    "for i in skeleton_ids:\n",
    "    if (i >2000)&(i<3000):\n",
    "        if i in Blue_cone_ids:\n",
    "            cell_ids.append([i,i,81])\n",
    "        elif i in cones_outside_cbc9:\n",
    "            cell_ids.append([i,i,82])\n",
    "        elif i in cut_off_cones:\n",
    "            cell_ids.append([i,i,83])\n",
    "        else:\n",
    "            cell_ids.append([i,i,80])\n",
    "    elif (i>3000):\n",
    "        cell_ids.append([i,i,84])\n",
    "cell_ids=np.array(cell_ids).astype(int)\n",
    "\n",
    "cell_ids_alt=cell_ids.copy()\n",
    "cell_ids_alt[np.in1d(cell_ids_alt[:,0],Blue_cone_ids_alt),2]=81"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Determine image size\n",
    "# total_min=np.ones(2)*500\n",
    "# total_max=np.zeros(2)\n",
    "# for i in range(cell_ids.shape[0]):\n",
    "#     if cell_ids[i,2]>=80:\n",
    "#         try:\n",
    "#             data=((np.loadtxt('data/cone_projections/cell'+str(cell_ids[i,0]).zfill(4)+'_rot_x_complete.gz'))).astype(int)\n",
    "#             total_min=np.minimum(total_min,np.min(data,axis=0))\n",
    "#             total_max=np.maximum(total_max,np.max(data,axis=0))\n",
    "#         except FileNotFoundError:\n",
    "#             continue\n",
    "# im_size=((total_max-total_min)/np.array([3,2])).astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "im_size=np.array([2355,1518])\n",
    "total_min=np.array([ -19, -213])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "im1=Image.new('CMYK',tuple(im_size+1),(0,0,0,0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Draw all cones\n",
    "for i in range(cell_ids.shape[0]):\n",
    "    if (cell_ids[i,2]>=80) and (cell_ids[i,2]<84):\n",
    "        try:\n",
    "            data=((np.loadtxt('data/cone_projections/cell'+str(cell_ids[i,0]).zfill(4)+'_rot_x_complete.gz')-total_min)/np.array([3,2])).astype(int)\n",
    "            draw = ImageDraw.Draw(im1)\n",
    "            draw.point(data.flatten().tolist(),fill=(0,0,0,80))\n",
    "            del draw\n",
    "        except FileNotFoundError:\n",
    "            continue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Draw green cones\n",
    "for i in range(cell_ids.shape[0]):\n",
    "    if cell_ids[i,2]==80:\n",
    "        try:\n",
    "            data=((np.loadtxt('data/cone_projections/cell'+str(cell_ids[i,0]).zfill(4)+'_rot_x_complete.gz')-total_min)/np.array([3,2])).astype(int)\n",
    "            draw = ImageDraw.Draw(im1)\n",
    "            draw.point(data.flatten().tolist(),fill=(194,0,255,0))\n",
    "            del draw\n",
    "        except FileNotFoundError:\n",
    "            continue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Draw blue cones\n",
    "for i in range(cell_ids.shape[0]):\n",
    "    if cell_ids[i,2]==81:\n",
    "        try:\n",
    "            data=((np.loadtxt('/data/cone_projections/cell'+str(cell_ids[i,0]).zfill(4)+'_rot_x_complete.gz')-total_min)/np.array([3,2])).astype(int)\n",
    "            draw = ImageDraw.Draw(im1)\n",
    "            draw.point(data.flatten().tolist(),fill=(232,181,0,0))\n",
    "            del draw\n",
    "        except FileNotFoundError:\n",
    "            continue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Draw cut off cones\n",
    "for i in range(cell_ids.shape[0]):\n",
    "    if cell_ids[i,2]==83:\n",
    "        try:\n",
    "            data=((np.loadtxt('data/cone_projections/cell'+str(cell_ids[i,0]).zfill(4)+'_rot_x_complete.gz')-total_min)/np.array([3,2])).astype(int)\n",
    "            draw = ImageDraw.Draw(im1)\n",
    "            draw.point(data.flatten().tolist(),fill=(0,0,0,255))\n",
    "            del draw\n",
    "        except FileNotFoundError:\n",
    "            continue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Draw cones not cut off\n",
    "for i in range(cell_ids.shape[0]):\n",
    "    if (cell_ids[i,2]>=80) and (cell_ids[i,2]<83):\n",
    "        try:\n",
    "            data=((np.loadtxt('../data/cone_projections/cell'+str(cell_ids[i,0]).zfill(4)+'_rot_x_complete.gz')-total_min)/np.array([3,2])).astype(int)\n",
    "            draw = ImageDraw.Draw(im1)\n",
    "            draw.point(data.flatten().tolist(),fill=(0,0,0,80))\n",
    "            del draw\n",
    "        except FileNotFoundError:\n",
    "            continue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "CBC9_colors=[(255,0,0,0),(255,255,0,0),(0,255,0,0),(0,255,255,0),(255,0,255,0),(0,0,255,0)]\n",
    "CBC9=[690,692,693,694,695,696]\n",
    "\n",
    "XBC_colors=[(32,255,189,10)]\n",
    "XBC_BC=[607]\n",
    "XBC_cones=[2106]\n",
    "\n",
    "dendritic_field_size_colors=[(32,255,189,10),(40,170,255,12)]\n",
    "dendritic_field_size_BC=[687,633]#[685,633]\n",
    "dendritic_field_cones=[2064,2070,2069,2178,2173,2097,2067,2184,2098,2079,2109,2093,2103]#[2021,2099,2102,2028,2103,2148,2147,2079,2109,2093,2103]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "soma_position_CBC9=[]\n",
    "for k in range(len(CBC9)):\n",
    "    i=cell_ids[cell_ids[:,0]==CBC9[k],1]\n",
    "    soma_position_CBC9.append(soma_pos[soma_line_ids[0,np.where(soma_internal_ids==i)[1][0]]-1,:])\n",
    "soma_position_CBC9=np.array(soma_position_CBC9)\n",
    "soma_position_CBC9[:,1:]-=total_min"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Draw type 9 BC\n",
    "draw = ImageDraw.Draw(im1)\n",
    "for k in range(len(CBC9)):\n",
    "    for cell in np.where(skeleton_ids==CBC9[k])[0]:\n",
    "        if skeletons[cell].item()[0].shape[1]==4:\n",
    "            if skeletons[cell].item()[0].shape[0]<5:\n",
    "                continue\n",
    "            nodes=np.copy(skeletons[cell].item()[0])\n",
    "            edges=np.copy(skeletons[cell].item()[2])\n",
    "        elif skeletons[cell].item()[2].shape[1]==4:\n",
    "            nodes=np.copy(skeletons[cell].item()[2])\n",
    "            edges=np.copy(skeletons[cell].item()[4])\n",
    "        else:\n",
    "            print('Index problem at',i)\n",
    "            break\n",
    "        nodes[:,:3]=nodes[:,:3]/[16.5,16.5,25]\n",
    "        \n",
    "        nodes[:,:3]=np.dot(M,nodes[:,:3].T).T\n",
    "        for i in range(edges.shape[0]):\n",
    "            if (nodes[edges[i,0]-1,0]<soma_position_CBC9[k,0]) and (nodes[edges[i,1]-1,0]<soma_position_CBC9[k,0]):\n",
    "                draw.line([tuple(nodes[edges[i,0]-1,1:3]/[3,2]-total_min/[3,2]),tuple(nodes[edges[i,1]-1,1:3]/[3,2]-total_min/[3,2])],fill=CBC9_colors[k],width=5)\n",
    "del draw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#draw CBC9 somata\n",
    "draw=ImageDraw.Draw(im1)\n",
    "r=30\n",
    "for k in range(len(CBC9)):\n",
    "    draw.ellipse((soma_position_CBC9[k,1]/3-r, soma_position_CBC9[k,2]/2-r, soma_position_CBC9[k,1]/3+r, soma_position_CBC9[k,2]/2+r), fill=CBC9_colors[k])\n",
    "del draw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Draw 10µm scale bar\n",
    "draw = ImageDraw.Draw(im1)\n",
    "draw.rectangle([2100,1450,2300,1470],fill=(0,0,0,255))\n",
    "del draw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Draw image of dendritic field size (BC)\n",
    "for cone in dendritic_field_cones:\n",
    "    try:\n",
    "        data=((np.loadtxt('data/cone_projections/cell'+str(cone).zfill(4)+'_rot_x_complete.gz')-total_min)/np.array([3,2])).astype(int)\n",
    "        draw = ImageDraw.Draw(im1)\n",
    "        draw.point(data.flatten().tolist(),fill=(100,0,0,80))\n",
    "        del draw\n",
    "    except FileNotFoundError:\n",
    "        continue\n",
    "\n",
    "draw = ImageDraw.Draw(im1)\n",
    "for k in range(len(dendritic_field_size_BC)):\n",
    "    for cell in np.where(skeleton_ids==dendritic_field_size_BC[k])[0]:\n",
    "        if skeletons[cell].item()[0].shape[1]==4:\n",
    "            if skeletons[cell].item()[0].shape[0]<5:\n",
    "                continue\n",
    "            nodes=np.copy(skeletons[cell].item()[0])\n",
    "            edges=np.copy(skeletons[cell].item()[2])\n",
    "        elif skeletons[cell].item()[2].shape[1]==4:\n",
    "            nodes=np.copy(skeletons[cell].item()[2])\n",
    "            edges=np.copy(skeletons[cell].item()[4])\n",
    "        else:\n",
    "            print('Index problem at',i)\n",
    "            break\n",
    "        nodes[:,:3]=nodes[:,:3]/[16.5,16.5,25]\n",
    "        nodes[:,:3]=np.dot(M,nodes[:,:3].T).T\n",
    "        for i in range(edges.shape[0]):\n",
    "            if (nodes[edges[i,0]-1,0]<4000) and (nodes[edges[i,1]-1,0]<4000):\n",
    "                draw.line([tuple(nodes[edges[i,0]-1,1:3]/[3,2]-total_min/[3,2]),tuple(nodes[edges[i,1]-1,1:3]/[3,2]-total_min/[3,2])],fill=dendritic_field_size_colors[k],width=5)\n",
    "del draw\n",
    "im2=im1.crop((360,418,1660,1518))\n",
    "draw = ImageDraw.Draw(im2)#5µm scale bar\n",
    "draw.rectangle([1050,1030,1250,1050],fill=(0,0,0,255))#10µm scale bar\n",
    "del draw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Draw image of XBC\n",
    "for cone in XBC_cones:\n",
    "    try:\n",
    "        data=((np.loadtxt('data/cone_projections/cell'+str(cone).zfill(4)+'_rot_x_complete.gz')-total_min)/np.array([3,2])).astype(int)\n",
    "        draw = ImageDraw.Draw(im1)\n",
    "        draw.point(data.flatten().tolist(),fill=(100,0,0,80))\n",
    "        del draw\n",
    "    except FileNotFoundError:\n",
    "        continue\n",
    "        \n",
    "for cone in [2006,2036]:\n",
    "    try:\n",
    "        data=((np.loadtxt('data/cone_projections/cell'+str(cone).zfill(4)+'_rot_x_complete.gz')-total_min)/np.array([3,2])).astype(int)\n",
    "        draw = ImageDraw.Draw(im1)\n",
    "        draw.point(data.flatten().tolist(),fill=(0,0,0,150))\n",
    "        del draw\n",
    "    except FileNotFoundError:\n",
    "        continue\n",
    "        \n",
    "draw = ImageDraw.Draw(im1)\n",
    "for k in range(len(XBC_BC)):\n",
    "    for cell in np.where(skeleton_ids==XBC_BC[k])[0]:\n",
    "        if skeletons[cell].item()[0].shape[1]==4:\n",
    "            if skeletons[cell].item()[0].shape[0]<5:\n",
    "                continue\n",
    "            nodes=np.copy(skeletons[cell].item()[0])\n",
    "            edges=np.copy(skeletons[cell].item()[2])\n",
    "        elif skeletons[cell].item()[2].shape[1]==4:\n",
    "            nodes=np.copy(skeletons[cell].item()[2])\n",
    "            edges=np.copy(skeletons[cell].item()[4])\n",
    "        else:\n",
    "            print('Index problem at',i)\n",
    "            break\n",
    "        nodes[:,:3]=nodes[:,:3]/[16.5,16.5,25]\n",
    "        nodes[:,:3]=np.dot(M,nodes[:,:3].T).T\n",
    "        for i in range(edges.shape[0]):\n",
    "            if (nodes[edges[i,0]-1,0]<4000) and (nodes[edges[i,1]-1,0]<4000):\n",
    "                draw.line([tuple(nodes[edges[i,0]-1,1:3]/[3,2]-total_min/[3,2]),tuple(nodes[edges[i,1]-1,1:3]/[3,2]-total_min/[3,2])],fill=XBC_colors[k],width=5)\n",
    "del draw\n",
    "im2=im1.crop((1250,300,1900,1000))\n",
    "draw = ImageDraw.Draw(im2)\n",
    "draw.rectangle([50,640,150,650],fill=(0,0,0,255))#5µm scale bar\n",
    "del draw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(20,10))\n",
    "plt.imshow(im1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# im1.save('figures/cone_overview.tiff')"
   ]
  },
  {
   "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
}