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 [ ]: