Cone projection images

This notebook contains the code to draw all images showing projections of the cone volumes along the optical axis such as figures 1C&D


In [ ]:
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from scipy.io import loadmat
from scipy.spatial import ConvexHull
from PIL import Image
from PIL import ImageDraw
from PIL import ImageFont
from shapely.geometry import Point
from shapely.geometry import Polygon
import pandas as pd
%matplotlib inline
matplotlib.rcParams.update({'font.size': 14})

In [ ]:
matplotlib.rcParams['font.family'] = 'sans-serif'
matplotlib.rcParams['font.sans-serif'] = ['Arial']

In [ ]:
skeleton_data=loadmat('data/skeletons_OPL_final.mat')
skeleton_ids=skeleton_data['kn_allSkeletons_e2006_IDs'].flatten()
skeletons=skeleton_data['kn_allSkeletons_e2006'].flatten()

In [ ]:
M=np.loadtxt('data/coordinate_rotation')

In [ ]:
soma_data=loadmat('../data/Soma_positions.mat')
soma_pos=soma_data['kn_e2006_ALLSKELETONS_FINAL2012_allSomata']
soma_internal_ids=soma_data['kn_e2006_ALLSKELETONS_FINAL2012_cellIDs']
soma_line_ids=soma_data['kn_e2006_ALLSKELETONS_FINAL2012_cellIDs_pure_forSomata']
soma_pos=np.dot(M,soma_pos[:,:3].T).T

In [ ]:
cut_off_cones=np.array([73,72,71,74,201,202,211,203,88,204,142,34,4,30,199])+2000
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
Blue_cone_ids=np.array([2006,2007,2009,2024,2028,2029])
Blue_cone_ids_alt=np.array([2007,2006,2010,2029,2061,2009,2031,2014,2024,2023,2100,2028,2017,2090])

In [ ]:
BC_ids=np.loadtxt('data/BC_IDs_new')

cell_ids=[]
for i in range(BC_ids.shape[0]):
    if BC_ids[i,4] in range(58,72):
        cell_ids.append([BC_ids[i,0],BC_ids[i,1],BC_ids[i,4]])
for i in skeleton_ids:
    if (i >2000)&(i<3000):
        if i in Blue_cone_ids:
            cell_ids.append([i,i,81])
        elif i in cones_outside_cbc9:
            cell_ids.append([i,i,82])
        elif i in cut_off_cones:
            cell_ids.append([i,i,83])
        else:
            cell_ids.append([i,i,80])
    elif (i>3000):
        cell_ids.append([i,i,84])
cell_ids=np.array(cell_ids).astype(int)

cell_ids_alt=cell_ids.copy()
cell_ids_alt[np.in1d(cell_ids_alt[:,0],Blue_cone_ids_alt),2]=81

In [ ]:
# Determine image size
# total_min=np.ones(2)*500
# total_max=np.zeros(2)
# for i in range(cell_ids.shape[0]):
#     if cell_ids[i,2]>=80:
#         try:
#             data=((np.loadtxt('data/cone_projections/cell'+str(cell_ids[i,0]).zfill(4)+'_rot_x_complete.gz'))).astype(int)
#             total_min=np.minimum(total_min,np.min(data,axis=0))
#             total_max=np.maximum(total_max,np.max(data,axis=0))
#         except FileNotFoundError:
#             continue
# im_size=((total_max-total_min)/np.array([3,2])).astype(int)

In [ ]:
im_size=np.array([2355,1518])
total_min=np.array([ -19, -213])

In [ ]:
im1=Image.new('CMYK',tuple(im_size+1),(0,0,0,0))

In [ ]:
#Draw all cones
for i in range(cell_ids.shape[0]):
    if (cell_ids[i,2]>=80) and (cell_ids[i,2]<84):
        try:
            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)
            draw = ImageDraw.Draw(im1)
            draw.point(data.flatten().tolist(),fill=(0,0,0,80))
            del draw
        except FileNotFoundError:
            continue

In [ ]:
#Draw green cones
for i in range(cell_ids.shape[0]):
    if cell_ids[i,2]==80:
        try:
            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)
            draw = ImageDraw.Draw(im1)
            draw.point(data.flatten().tolist(),fill=(194,0,255,0))
            del draw
        except FileNotFoundError:
            continue

In [ ]:
#Draw blue cones
for i in range(cell_ids.shape[0]):
    if cell_ids[i,2]==81:
        try:
            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)
            draw = ImageDraw.Draw(im1)
            draw.point(data.flatten().tolist(),fill=(232,181,0,0))
            del draw
        except FileNotFoundError:
            continue

In [ ]:
#Draw cut off cones
for i in range(cell_ids.shape[0]):
    if cell_ids[i,2]==83:
        try:
            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)
            draw = ImageDraw.Draw(im1)
            draw.point(data.flatten().tolist(),fill=(0,0,0,255))
            del draw
        except FileNotFoundError:
            continue

In [ ]:
#Draw cones not cut off
for i in range(cell_ids.shape[0]):
    if (cell_ids[i,2]>=80) and (cell_ids[i,2]<83):
        try:
            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)
            draw = ImageDraw.Draw(im1)
            draw.point(data.flatten().tolist(),fill=(0,0,0,80))
            del draw
        except FileNotFoundError:
            continue

In [ ]:
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)]
CBC9=[690,692,693,694,695,696]

XBC_colors=[(32,255,189,10)]
XBC_BC=[607]
XBC_cones=[2106]

dendritic_field_size_colors=[(32,255,189,10),(40,170,255,12)]
dendritic_field_size_BC=[687,633]#[685,633]
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]

In [ ]:
soma_position_CBC9=[]
for k in range(len(CBC9)):
    i=cell_ids[cell_ids[:,0]==CBC9[k],1]
    soma_position_CBC9.append(soma_pos[soma_line_ids[0,np.where(soma_internal_ids==i)[1][0]]-1,:])
soma_position_CBC9=np.array(soma_position_CBC9)
soma_position_CBC9[:,1:]-=total_min

In [ ]:
#Draw type 9 BC
draw = ImageDraw.Draw(im1)
for k in range(len(CBC9)):
    for cell in np.where(skeleton_ids==CBC9[k])[0]:
        if skeletons[cell].item()[0].shape[1]==4:
            if skeletons[cell].item()[0].shape[0]<5:
                continue
            nodes=np.copy(skeletons[cell].item()[0])
            edges=np.copy(skeletons[cell].item()[2])
        elif skeletons[cell].item()[2].shape[1]==4:
            nodes=np.copy(skeletons[cell].item()[2])
            edges=np.copy(skeletons[cell].item()[4])
        else:
            print('Index problem at',i)
            break
        nodes[:,:3]=nodes[:,:3]/[16.5,16.5,25]
        
        nodes[:,:3]=np.dot(M,nodes[:,:3].T).T
        for i in range(edges.shape[0]):
            if (nodes[edges[i,0]-1,0]<soma_position_CBC9[k,0]) and (nodes[edges[i,1]-1,0]<soma_position_CBC9[k,0]):
                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)
del draw

In [ ]:
#draw CBC9 somata
draw=ImageDraw.Draw(im1)
r=30
for k in range(len(CBC9)):
    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])
del draw

In [ ]:
#Draw 10µm scale bar
draw = ImageDraw.Draw(im1)
draw.rectangle([2100,1450,2300,1470],fill=(0,0,0,255))
del draw

In [ ]:
#Draw image of dendritic field size (BC)
for cone in dendritic_field_cones:
    try:
        data=((np.loadtxt('data/cone_projections/cell'+str(cone).zfill(4)+'_rot_x_complete.gz')-total_min)/np.array([3,2])).astype(int)
        draw = ImageDraw.Draw(im1)
        draw.point(data.flatten().tolist(),fill=(100,0,0,80))
        del draw
    except FileNotFoundError:
        continue

draw = ImageDraw.Draw(im1)
for k in range(len(dendritic_field_size_BC)):
    for cell in np.where(skeleton_ids==dendritic_field_size_BC[k])[0]:
        if skeletons[cell].item()[0].shape[1]==4:
            if skeletons[cell].item()[0].shape[0]<5:
                continue
            nodes=np.copy(skeletons[cell].item()[0])
            edges=np.copy(skeletons[cell].item()[2])
        elif skeletons[cell].item()[2].shape[1]==4:
            nodes=np.copy(skeletons[cell].item()[2])
            edges=np.copy(skeletons[cell].item()[4])
        else:
            print('Index problem at',i)
            break
        nodes[:,:3]=nodes[:,:3]/[16.5,16.5,25]
        nodes[:,:3]=np.dot(M,nodes[:,:3].T).T
        for i in range(edges.shape[0]):
            if (nodes[edges[i,0]-1,0]<4000) and (nodes[edges[i,1]-1,0]<4000):
                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)
del draw
im2=im1.crop((360,418,1660,1518))
draw = ImageDraw.Draw(im2)#5µm scale bar
draw.rectangle([1050,1030,1250,1050],fill=(0,0,0,255))#10µm scale bar
del draw

In [ ]:
#Draw image of XBC
for cone in XBC_cones:
    try:
        data=((np.loadtxt('data/cone_projections/cell'+str(cone).zfill(4)+'_rot_x_complete.gz')-total_min)/np.array([3,2])).astype(int)
        draw = ImageDraw.Draw(im1)
        draw.point(data.flatten().tolist(),fill=(100,0,0,80))
        del draw
    except FileNotFoundError:
        continue
        
for cone in [2006,2036]:
    try:
        data=((np.loadtxt('data/cone_projections/cell'+str(cone).zfill(4)+'_rot_x_complete.gz')-total_min)/np.array([3,2])).astype(int)
        draw = ImageDraw.Draw(im1)
        draw.point(data.flatten().tolist(),fill=(0,0,0,150))
        del draw
    except FileNotFoundError:
        continue
        
draw = ImageDraw.Draw(im1)
for k in range(len(XBC_BC)):
    for cell in np.where(skeleton_ids==XBC_BC[k])[0]:
        if skeletons[cell].item()[0].shape[1]==4:
            if skeletons[cell].item()[0].shape[0]<5:
                continue
            nodes=np.copy(skeletons[cell].item()[0])
            edges=np.copy(skeletons[cell].item()[2])
        elif skeletons[cell].item()[2].shape[1]==4:
            nodes=np.copy(skeletons[cell].item()[2])
            edges=np.copy(skeletons[cell].item()[4])
        else:
            print('Index problem at',i)
            break
        nodes[:,:3]=nodes[:,:3]/[16.5,16.5,25]
        nodes[:,:3]=np.dot(M,nodes[:,:3].T).T
        for i in range(edges.shape[0]):
            if (nodes[edges[i,0]-1,0]<4000) and (nodes[edges[i,1]-1,0]<4000):
                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)
del draw
im2=im1.crop((1250,300,1900,1000))
draw = ImageDraw.Draw(im2)
draw.rectangle([50,640,150,650],fill=(0,0,0,255))#5µm scale bar
del draw

In [ ]:
plt.figure(figsize=(20,10))
plt.imshow(im1)
plt.show()

In [ ]:
# im1.save('figures/cone_overview.tiff')

In [ ]: