"""
Basic visualization of neurite morphologies using matplotlib.
Usage is restricted to morphologies in the sWC format with the three-point soma `standard <http://neuromorpho.org/neuroMorpho/SomaFormat.html>`_
B. Torben-Nielsen
"""
import sys,time
sys.setrecursionlimit(10000)
import numpy as np
import math
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import matplotlib.cm as cm
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.gridspec as gridspec
import btmorph
from btmorph import config
from numpy import mean,cov,double,cumsum,dot,linalg,array,rank
from pylab import plot,subplot,axis,stem,show,figure
""" internal constants required for the dendrogram generation """
H_SPACE = 20
V_SPACE = 0
C = 'k'
max_width = 0
max_height = 0
def _plot_3D_figure(SWC,my_color_list = ['r','g','b','c','m','y','r--','b--','g--']):
"""
Parameters
-----------
SWC: dict
"""
for index in SWC.keys() : # not ordered but that has little importance here
# draw a line segment from parent to current point
current_SWC = SWC[index]
#print 'index: ', index, ' -> ', current_SWC
c_x = current_SWC[0]
c_y = current_SWC[1]
c_z = current_SWC[2]
c_r = current_SWC[3]
parent_index = current_SWC[4]
if(index <= 3) :
#print 'do not draw the soma and its CNG, !!! 2 !!! point descriptions'
pass
else :
parent_SWC = SWC[parent_index]
p_x = parent_SWC[0]
p_y = parent_SWC[1]
p_z = parent_SWC[2]
p_r = parent_SWC[3]
# print 'index:', index, ', len(cs)=', len(cs)
pl = plt.plot([p_x,c_x],[p_y,c_y],[p_z,c_z],my_color_list[current_SWC[5]-1],linewidth=c_r)
[docs]def true_2D_projections(file_name='P20-DEV139.CNG.swc',align=True,outN=None,bar=None,depth="Z"):
"""
Three 2D projections
Parameters
-----------
file_name : string
File name of the SWC file to plots
depth : string
Set which axis represents "depth". In experimental work, the \
Z-axis is depth (as in my PPNeurMorphC) but in NeuroMorpho the \
Y-axis is the depth. (Depth is used to denote the axis from deep to superficial)
align : boolean
Translate the figure so that the soma is on the origin [0,0,0].
outN : string
File name of the output file. Extension of this file sets the file type
bar : list of int or real
Three values to set the thicks and marks on the plot in X,Y and Z-dimension
depth : string
Set the axis representing the depth (= axis from superficial to deep). In most SWC files \
this is 'Y'. The other option is 'Z', that is more consistent with the usual Cartesian coordinate systems
"""
my_color_list = ['r','g','b','c','m','y','k','g','DarkGray']
# read the SWC into a dictionary: key=index, value=(x,y,z,d,parent)
x = open(file_name,'r')
soma = None
SWC = {}
for line in x :
if(not line.startswith('#')) :
splits = line.split()
index = int(splits[0])
if index == 1:
soma_x = float(splits[2])
soma_y = float(splits[3])
soma_z = float(splits[4])
n_type = int(splits[1])
if not align:
x = float(splits[2])
y = float(splits[3])
z = float(splits[4])
else:
x = float(splits[2])-soma_x
y = float(splits[3])-soma_y
z = float(splits[4])-soma_z
r = float(splits[5])
parent = int(splits[-1])
SWC[index] = (x,y,z,r,parent,n_type)
fig = plt.figure()
if depth == "Y":
ax = plt.subplot2grid((2,2), (0, 0))
for index in SWC.keys():
if index < 3:
continue
C = SWC[index]
P = SWC[C[4]] # C[4] is parent index
pl = plt.plot([P[0],C[0]],[P[1],C[1]],my_color_list[C[5]-1],linewidth=C[3])
plt.xlabel("X")
plt.ylabel("Y")
plt.xticks([0,bar[0]])
plt.yticks([0,bar[1]])
ax2 = plt.subplot2grid((2,2), (0, 1))
for index in SWC.keys():
if index < 3:
continue
C = SWC[index]
P = SWC[C[4]] # C[4] is parent index
pl2 = plt.plot([P[2],C[2]],[P[1],C[1]],my_color_list[C[5]-1],linewidth=C[3])
plt.xlabel("Z")
plt.ylabel("Y")
plt.xticks([0,bar[2]])
plt.yticks([0,bar[1]])
ax3 = plt.subplot2grid((2,2), (1, 1))
for index in SWC.keys():
if index < 3:
continue
C = SWC[index]
P = SWC[C[4]] # C[4] is parent index
pl3 = plt.plot([P[0],C[0]],[P[2],C[2]],my_color_list[C[5]-1],linewidth=C[3])
plt.xlabel("X")
plt.ylabel("Z")
plt.xticks([0,bar[0]])
plt.yticks([0,bar[2]])
else: # default for my code: Z is depth
ax = plt.subplot2grid((2,2), (0, 0))
for index in SWC.keys():
if index < 3:
continue
C = SWC[index]
P = SWC[C[4]] # C[4] is parent index
pl = plt.plot([P[0],C[0]],[P[2],C[2]],my_color_list[C[5]-1],linewidth=C[3])
plt.xlabel("X")
plt.ylabel("Z")
plt.xticks([0,bar[0]])
plt.yticks([0,bar[2]])
ax2 = plt.subplot2grid((2,2), (0, 1))
for index in SWC.keys():
if index < 3:
continue
C = SWC[index]
P = SWC[C[4]] # C[4] is parent index
pl2 = plt.plot([P[1],C[1]],[P[2],C[2]],my_color_list[C[5]-1],linewidth=C[3])
plt.xlabel("Y")
plt.ylabel("Z")
plt.xticks([0,bar[1]])
plt.yticks([0,bar[2]])
ax3 = plt.subplot2grid((2,2), (1, 1))
for index in SWC.keys():
if index < 3:
continue
C = SWC[index]
P = SWC[C[4]] # C[4] is parent index
pl3 = plt.plot([P[0],C[0]],[P[1],C[1]],my_color_list[C[5]-1],linewidth=C[3])
plt.xlabel("X")
plt.ylabel("Y")
plt.xticks([0,bar[0]])
plt.yticks([0,bar[1]])
plt.tight_layout()
if not outN == None:
plt.savefig(outN)
[docs]def plot_2D_SWC(file_name='P20-DEV139.CNG.swc',cs=None,\
synapses=None,syn_cs='ro',\
outN=None,align=True,offset=None,\
show_axis=False,depth="Y",filter=range(10),\
show_radius=True,bar_L=None,bar=[50,50,50],\
new_fig=True,color_scheme="default"):
"""
2D matplotlib plot of a neuronal moprhology. Projection can be in XY and XZ.
The SWC has to be formatted with a "three point soma".
Colors can be provided
Parameters
-----------
file_name : string
File name of the SWC file to plots
cs : list of float
Raw values that will be plotted on the morphology according to a colormap
synapses : list of int
Indices of the compartments where synapses (= small circles) should be drawn
syn_c : string
Color of the synapses ('r','g', 'b', ...). String follows matplotlib conventions. You can \
include the marker style as well. Default `syn_c='ro'`
outN : string
File name of the output file. Extension of this file sets the file type
align : boolean
Translate the figure so that the soma is on the origin [0,0,0].
offset : list on float
List of length 3 with X,Y and Z shift of the morphology to be plotted. Not to be used in conjunction with the "align" option
show_axis : boolean
whether or not to draw the axis
filter : list
List of integers indicating the SWC types to be included (1:soma, 2:axon, 3:basal,4:apical,...). By default, all SWC types are included
show_radius : boolean
Default "True" plots the diameter; "False" will render a wire plot.
bar_L : float
Add a scale bar to the plot. *Currently only works with align=True*
bar : list of real
When the axis are shown (`show_axis=True`), ticks will be plotted according to this list.\
List contains 3 values for X, Y and Z ticks. Default [50,50,50]
depth : string
Default "Y" means that X represents the superficial to deep axis. \
Otherwise, use "Z" to conform the mathematical standard of having the Z axis.
new_fig : boolean
True if matplotlib has to plot in a new figure. False, if otherwise.
color_scheme : string
Set the color scheme used for background and neurite types. Default `default`. Other option `neuromorpho`
Notes
-----
If the soma is not located at [0,0,0], the scale bar (`bar_L`) and the ticks (`bar`) might not work \
as expected
"""
#print "scheme: ", config.c_scheme_nm
plot_radius = 0.5
#my_color_list = ['r','g','b','c','m','y','k','g','DarkGray']
if color_scheme == 'default':
my_color_list = config.c_scheme_default['neurite']
elif color_scheme == 'neuromorpho':
my_color_list = config.c_scheme_nm['neurite']
print 'my_color_list: ', my_color_list
# resolve some potentially conflicting arguments
if not offset == None:
align = False
# read the SWC into a dictionary: key=index, value=(x,y,z,d,parent)
x = open(file_name,'r')
soma = None
SWC = {}
for line in x :
if(not line.startswith('#')) :
splits = line.split()
index = int(splits[0])
if index == 1:
if offset == None:
soma_x = float(splits[2])
soma_y = float(splits[3])
soma_z = float(splits[4])
else:
soma_x = float(splits[2]) + offset[0]
soma_y = float(splits[3]) + offset[1]
soma_z = float(splits[4]) + offset[2]
n_type = int(splits[1])
if not align:
if offset == None:
x = float(splits[2])
y = float(splits[3])
z = float(splits[4])
else:
x = float(splits[2]) + offset[0]
y = float(splits[3]) + offset[1]
z = float(splits[4]) + offset[2]
else:
x = float(splits[2])-soma_x
y = float(splits[3])-soma_y
z = float(splits[4])-soma_z
r = float(splits[5])
parent = int(splits[-1])
if n_type in filter:
SWC[index] = (x,y,z,r,parent,n_type)
# setting up for a colormap
if cs == None :
pass
else :
norm = colors.normalize(np.min(cs),np.max(cs))
Z = [[0,0],[0,0]]
levels=range(int(np.min(cs)),int(np.max(cs)),100)
levels = np.linspace(np.min(cs),np.max(cs),100)
CS3 = plt.contourf(Z,levels,cmap=cm.jet)
plt.clf()
if new_fig:
fig = plt.figure()
min_depth=100
if depth == "Y":
#ax = plt.subplot2grid((2,2), (0, 0))
for index in SWC.keys():
if index < 3:
continue
C = SWC[index]
P = SWC[C[4]] # C[4] is parent index
line_width = C[3] if show_radius else plot_radius
if cs == None:
pl = plt.plot([P[0],C[0]],[P[1],C[1]],my_color_list[C[5]-1],linewidth=line_width)
else:
#c=cm.jet(norm(cs[index]))
pl = plt.plot([P[0],C[0]],[P[1],C[1]],c=cm.jet(norm(cs[index])),linewidth=line_width)
# add the synapses
if not synapses == None :
if index in synapses :
plt.plot([C[0]],C[1],syn_cs)
min_depth = C[1] if C[1] < min_depth else min_depth
# TODO: insert synapses here
plt.xlabel("X")
plt.ylabel("Y")
plt.xticks([0,bar[0]])
plt.yticks([0,bar[1]])
else: # default for my code: Z is depth
#ax = plt.subplot2grid((2,2), (0, 0))
for index in SWC.keys():
if index < 3:
continue
C = SWC[index]
P = SWC[C[4]] # C[4] is parent index
line_width = C[3] if show_radius else plot_radius
pl = plt.plot([P[0],C[0]],[P[2],C[2]],my_color_list[C[5]-1],linewidth=line_width)
# add the synapses
if(synapses != None) :
if index in synapses :
plt.plot([C[0]],C[2],syn_cs)
min_depth = C[2] if C[2] < min_depth else min_depth
# TODO: insert synapses here
plt.xlabel("X")
plt.ylabel("Z")
plt.xticks([0,bar[0]])
plt.yticks([0,bar[2]])
plt.tight_layout()
plt.axis("equal")
# axes or not?
frame1 = plt.gca()
if not show_axis :
frame1.axes.get_xaxis().set_ticks([])
frame1.axes.get_yaxis().set_ticks([])
frame1.axes.get_xaxis().set_visible(False)
frame1.axes.get_yaxis().set_visible(False)
# scale bar or not?
if offset == None and not bar_L == None:
plt.plot([0,bar_L],[min_depth*1.1,min_depth*1.1],'k',linewidth=5)
else :
pass # when offset is used, the figure is to scale and no scalebar needed
# color bar? in case `cs` is used
if not cs ==None :
cb = plt.colorbar(CS3) # bit of a workaround, but it seems to work
ticks_f = np.linspace(np.min(cs)-1,np.max(cs)+1,5)
ticks_i = map(int,ticks_f)
print "ticks_i: ", ticks_i
cb.set_ticks(ticks_i)
# set the bg color
fig = plt.gcf()
ax = fig.gca()
if color_scheme == 'default':
ax.set_axis_bgcolor(config.c_scheme_default['bg'])
elif color_scheme == 'neuromorpho':
ax.set_axis_bgcolor(config.c_scheme_nm['bg'])
if not outN == None:
plt.savefig(outN)
[docs]def plot_3D_SWC(file_name='P20-DEV139.CNG.swc',cs=None,synapses=None,syn_cs=None,outN=None,offset=None,align=True,filter=range(10)) :
"""
3D matplotlib plot of a neuronal morphology. The SWC has to be formatted with a "three point soma".
Colors can be provided and synapse location marked
Parameters
-----------
file_name : string
File name of the SWC file to plots
cs : list of float
Raw values that will be plotted on the morphology according to a colormap
synapses : list of int
Indices of the compartments where synapses (= small circles) should be drawn
syn_cs : string
Color of the synapses ('r','g', 'b', ...)
outN : string
File name of the output file. Extension of this file sets the file type
offset : list on float
List of length 3 with X,Y and Z shift of the morphology to be plotted. Not to be used in conjunction with the "align" option
show_axis : boolean
whether or not to draw the axis
filter : list
List of integers indicating the SWC types to be included (1:soma, 2:axon, 3:basal,4:apical,...). By default, all SWC types are included
"""
my_color_list = ['r','g','b','c','m','y','r--','b--','g--']
# resolve some potentially conflicting arguments
if not offset == None:
align = False
if(cs == None) :
pass
else :
norm = colors.normalize(np.min(cs),np.max(cs))
Z = [[0,0],[0,0]]
levels=range(int(np.min(cs)),int(np.max(cs)),100)
levels = np.linspace(np.min(cs),np.max(cs),100)
CS3 = plt.contourf(Z,levels,cmap=cm.jet)
plt.clf()
# read the SWC into a dictionary: key=index, value=(x,y,z,d,parent)
x = open(file_name,'r')
SWC = {}
# for line in x :
# if(not line.startswith('#')) :
# splits = line.split()
# index = int(splits[0])
# n_type = int(splits[1])
# x = float(splits[2])
# y = float(splits[3])
# z = float(splits[4])
# r = float(splits[5])
# parent = int(splits[-1])
# SWC[index] = (x,y,z,r,parent,n_type)
for line in x :
if(not line.startswith('#')) :
splits = line.split()
index = int(splits[0])
if index == 1:
soma_x = float(splits[2])
soma_y = float(splits[3])
soma_z = float(splits[4])
n_type = int(splits[1])
if not offset == None :
x = float(splits[2])-offset[0]
y = float(splits[3])-offset[1]
z = float(splits[4])-offset[2]
elif not align:
x = float(splits[2])
y = float(splits[3])
z = float(splits[4])
else:
x = float(splits[2])-soma_x
y = float(splits[3])-soma_y
z = float(splits[4])-soma_z
r = float(splits[5])
parent = int(splits[-1])
if n_type in filter:
SWC[index] = (x,y,z,r,parent,n_type)
fig = plt.figure()
ax = fig.gca(projection='3d')
for index in SWC.keys() : # not ordered but that has little importance here
# draw a line segment from parent to current point
current_SWC = SWC[index]
#print 'index: ', index, ' -> ', current_SWC
c_x = current_SWC[0]
c_y = current_SWC[1]
c_z = current_SWC[2]
c_r = current_SWC[3]
parent_index = current_SWC[4]
if(index <= 3) :
print 'do not draw the soma and its CNG, !!! 2 !!! point descriptions'
else :
parent_SWC = SWC[parent_index]
p_x = parent_SWC[0]
p_y = parent_SWC[1]
p_z = parent_SWC[2]
p_r = parent_SWC[3]
# print 'index:', index, ', len(cs)=', len(cs)
if cs == None :
pl = plt.plot([p_x,c_x],[p_y,c_y],[p_z,c_z],my_color_list[current_SWC[5]-1],linewidth=c_r/2.0)
else :
try :
pl = plt.plot([p_x,c_x],[p_y,c_y], \
c=cm.jet(norm(cs[index])),linewidth=c_r)
except Exception :
print 'something going wrong here'
# pass# it's ok: it's the list size...
# add the synapses
if(synapses != None) :
if(index in synapses) :
if syn_cs :
plt.plot(c_x,c_y,syn_cs)
else :
plt.plot(c_x,c_y,'ro')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
if(cs != None) :
cb = plt.colorbar(CS3) # bit of a workaround, but it seems to work
ticks_f = np.linspace(np.min(cs)-1,np.max(cs)+1,5)
ticks_i = map(int,ticks_f)
cb.set_ticks(ticks_i)
if(outN != None) :
plt.savefig(outN)
return ax
[docs]def plot_dendrogram(file_name,transform='plain',shift=0,c='k',radius=True,rm=20000.0,ra=200,outN=None) :
"""
Generate a dendrogram from an SWC file. The SWC has to be formatted with a "three point soma"
Parameters
-----------
file_name : string
File name of the SWC file to plots
transform : string
Either 'plain' or 'lambda'. Plain means no transform while 'lambda' performs an elecotrtonic transform
shift : float
Offset in the x-direction
c : string
Color ('r','g', 'b', ...)
radius : boolean
Plot a wire (False) dendrogram or one with the thickness of the processes (True)
rm : float
Membrane resistance. Only needed when transform = 'lambda'
rm : float
Axial resistance. Only needed when transform = 'lambda'
outN : string
File name of the output file. Extension of this file sets the file type
"""
global C, RM, RA, max_width, max_height # n.a.s.t.y.
swc_tree = btmorph.STree2()
swc_tree = swc_tree.read_SWC_tree_from_file(file_name)
RM = rm
RA = ra
C = c
max_height = 0
max_width = 0
plt.clf()
print 'Going to build the dendrogram. This might take a while...'
ttt = time.time()
_expand_dendrogram(swc_tree.root,swc_tree,shift,0,radius=radius,transform=transform)
if(transform == 'plain') :
plt.ylabel('L (micron)')
elif(transform == 'lambda') :
plt.ylabel('L (lambda)')
print (time.time() - ttt), ' later the dendrogram was finished. '
print 'max_widht=%f, max_height=%f' % (max_width,max_height)
x_bound = (max_width / 2.0) + (0.1*max_width)
max_y_bound = max_height + 0.1*max_height
plt.axis([-1.0*x_bound,x_bound,-0.1*max_height,max_y_bound])
plt.plot([x_bound,x_bound],[0,100],'k', linewidth=5) # 250 for MN, 100 for granule
frame1 = plt.gca()
frame1.axes.get_xaxis().set_visible(False)
frame1.axes.get_yaxis().set_visible(False)
if(outN != None) :
plt.savefig(outN)
def _expand_dendrogram(cNode,swc_tree,off_x,off_y,radius,transform='plain') :
global max_width,max_height # middle name d.i.r.t.y.
'''
Gold old fashioned recursion... sys.setrecursionlimit()!
'''
place_holder_h = H_SPACE
max_degree = swc_tree.degree_of_node(cNode)
required_h_space = max_degree * place_holder_h
start_x = off_x-(required_h_space/2.0)
if(required_h_space > max_width) :
max_width = required_h_space
if swc_tree.is_root(cNode) :
print 'i am expanding the root'
cNode.children.remove(swc_tree.get_node_with_index(2))
cNode.children.remove(swc_tree.get_node_with_index(3))
for cChild in cNode.children :
l = _path_between(swc_tree,cChild,cNode,transform=transform)
r = cChild.content['p3d'].radius
cChild_degree = swc_tree.degree_of_node(cChild)
new_off_x = start_x + ( (cChild_degree/2.0)*place_holder_h )
new_off_y = off_y+(V_SPACE*2)+l
r = r if radius else 1
plt.vlines(new_off_x,off_y+V_SPACE,new_off_y,linewidth=r,colors=C)
if((off_y+(V_SPACE*2)+l) > max_height) :
max_height = off_y+(V_SPACE*2)+l
_expand_dendrogram(cChild,swc_tree,new_off_x,new_off_y,radius=radius,transform=transform)
start_x = start_x + (cChild_degree*place_holder_h)
plt.hlines(off_y+V_SPACE,off_x,new_off_x,colors=C)
def _path_between(swc_tree,deep,high,transform='plain') :
path = swc_tree.path_to_root(deep)
pl = 0
pNode = deep
for node in path[1:] :
pPos = pNode.content['p3d'].xyz
cPos = node.content['p3d'].xyz
pl = pl + np.sqrt(np.sum((cPos-pPos)**2))
#pl += np.sqrt( (pPos.x - cPos.x)**2 + (pPos.y - cPos.y)**2 + (pPos.z - cPos.z)**2 )
pNode = node
if(node == high) : break
if(transform == 'plain'):
return pl
elif(transform == 'lambda') :
DIAM = (deep.content['p3d'].radius*2.0 + high.content['p3d'].radius*2.0) /2.0 # naive...
c_lambda = np.sqrt(1e+4*(DIAM/4.0)*(RM/RA))
return pl / c_lambda
def _pca(A):
""" performs principal components analysis
(PCA) on the n-by-p data matrix A
Rows of A correspond to observations, columns to variables.
Returns :
coeff :
is a p-by-p matrix, each column containing coefficients
for one principal component.
score :
the principal component scores; that is, the representation
of A in the principal component space. Rows of SCORE
correspond to observations, columns to components.
latent :
a vector containing the eigenvalues
of the covariance matrix of A.
source: http://glowingpython.blogspot.jp/2011/07/principal-component-analysis-with-numpy.html
"""
# computing eigenvalues and eigenvectors of covariance matrix
M = (A-mean(A.T,axis=1)).T # subtract the mean (along columns)
[latent,coeff] = linalg.eig(cov(M)) # attention:not always sorted
score = dot(coeff.T,M) # projection of the data in the new space
return coeff,score,latent
[docs]def pca_project_tree(tree):
"""
Returns a tree which is a projection of the original tree on the plane of most variance
Parameters
----------
tree : :class:`btmorph.btstructs2.STree2`
A tree
Returns
--------
tree : :class:`btmorph.btstructs2.STree2`
New flattened tree
"""
nodes = tree.get_nodes()
N = len(nodes)
coords = map(lambda n: n.content['p3d'].xyz, nodes)
points = np.transpose(coords)
coeff, score, latent = _pca(points.T)
score[2,:] = [0]*N
newp = np.transpose(score)
# Move soma to origin
translate = score[:,0]
for i in range(0, N):
nodes[i].content['p3d'].xyz = newp[i] - translate
import time
fmt = '%Y_%b_%d_%H_%M_%S'
now = time.strftime(fmt)
tree.write_SWC_tree_to_file('tmpTree_3d_'+now+ '.swc')
tree = btmorph.STree2().read_SWC_tree_from_file('tmpTree_3d_'+now+ '.swc')
import os
os.remove('tmpTree_3d_'+now+ '.swc')
return tree