BICCN PI Meeting Demo

Estimating Neuron Curvature/Torsion Demo

This notebook demonstrates fitting splines, and computing curvature/torsion from a sample neuron trace in SWC format

[1]:
from pathlib import Path
from brainlit.utils.Neuron_trace import NeuronTrace
from brainlit.algorithms.trace_analysis.fit_spline import GeometricGraph
from brainlit.algorithms.trace_analysis import spline_fxns
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
from scipy.interpolate import splev
#%matplotlib tk
/opt/buildhome/python3.7/lib/python3.7/site-packages/python_jsonschema_objects/__init__.py:53: UserWarning: Schema version http://json-schema.org/draft-04/schema not recognized. Some keywords and features may not be supported.
  self.schema["$schema"]

Supplemental function for visualization later

[2]:
def hist_equalize(array, bins):
    ra_histogram, bins = np.histogram(array, bins, density=True)
    cdf = ra_histogram.cumsum() # cumulative distribution function
    cdf = cdf / cdf[-1]
    # use linear interpolation of cdf to find new pixel values
    ra_equalized = np.interp(array, bins[:-1], cdf)

    return ra_equalized, cdf

Read SWC and fit splines

[3]:
brainlit_path = Path.cwd().parent.parent.parent
swc_path = Path.joinpath(brainlit_path,'data','data_octree','consensus-swcs','2018-08-01_G-002_consensus.swc')

nt = NeuronTrace(path=str(swc_path))
g = nt.get_graph()
df = nt.get_df()
neuron = GeometricGraph(df=df)
spline_tree = neuron.fit_spline_tree_invariant()
soma = np.array([df.x[0], df.y[0], df.z[0]])

View Splines

[4]:
def node_height(G, node):
    predecessors = list(G.predecessors(node))
    L = len(predecessors)
    assert L == 1 or L == 0
    if L == 0:
        return 0
    else:
        return 1 + node_height(G, predecessors[0])

fig = plt.figure(figsize=(12, 10), dpi=80)
ax = Axes3D(fig)
for edge in g.edges:
    n1 = g.nodes[edge[0]]
    n2 = g.nodes[edge[1]]
    ax.plot([n1['x'], n2['x']],[n1['y'], n2['y']],[n1['z'], n2['z']],'b', linewidth=0.5)


ax.set_axis_off()
ax.set_title("Piecewise Linear")
ax.view_init(-140, 60)

fig = plt.figure(figsize=(12, 10), dpi=80)
ax = Axes3D(fig)

for j, node in enumerate(spline_tree.nodes):
    spline = spline_tree.nodes[node]
    spline_height = node_height(spline_tree, node)
    tck, u_um = spline["spline"]
    y = splev(np.arange(u_um[0],u_um[-1], 0.1), tck)

    if spline_height == 0:
        c = "b"
        ax.scatter(y[0][0],y[1][0],y[2][0],'b', s=100)
    else:
        successors = spline_tree.successors(node)
        if len(list(successors)) == 0:
            c = "g"
        else:
            c = "r"

    ax.plot(y[0], y[1], y[2], c, linewidth=0.5)

ax.set_axis_off()
ax.set_title("Splines")
ax.view_init(-140, 60)
/opt/buildhome/python3.7/lib/python3.7/site-packages/ipykernel_launcher.py:11: MatplotlibDeprecationWarning: Axes3D(fig) adding itself to the figure is deprecated since 3.4. Pass the keyword argument auto_add_to_figure=False and use fig.add_axes(ax) to suppress this warning. The default value of auto_add_to_figure will change to False in mpl3.5 and True values will no longer work in 3.6.  This is consistent with other Axes classes.
  # This is added back by InteractiveShellApp.init_path()
/opt/buildhome/python3.7/lib/python3.7/site-packages/ipykernel_launcher.py:23: MatplotlibDeprecationWarning: Axes3D(fig) adding itself to the figure is deprecated since 3.4. Pass the keyword argument auto_add_to_figure=False and use fig.add_axes(ax) to suppress this warning. The default value of auto_add_to_figure will change to False in mpl3.5 and True values will no longer work in 3.6.  This is consistent with other Axes classes.
../../_images/notebooks_algorithms_biccn_demo_8_1.png
../../_images/notebooks_algorithms_biccn_demo_8_2.png

View one of the branches

[5]:
node = 7

path = spline_tree.nodes[node]["path"]
locs = np.zeros((len(path),3))
for p,point in enumerate(path):
    locs[p,:] = neuron.nodes[point]["loc"]

spline = spline_tree.nodes[node]["spline"]
u = spline[1]
u = np.arange(u[0], u[-1]+0.9, 1)
tck = spline[0]
pts = splev(u, tck)

fig = plt.figure(figsize=(12, 10), dpi=80)
ax = fig.add_subplot(111,projection="3d")
ax.plot(pts[0], pts[1], pts[2], 'red')
ax.w_xaxis.set_pane_color((0.23, 0.25, 0.209, 0.5))
ax.w_yaxis.set_pane_color((0.23, 0.25, 0.209, 0.1))
ax.w_zaxis.set_pane_color((0.23, 0.25, 0.209, 0.3))
ax.grid(False)
ax.set_xticks([])
ax.set_yticks([])
ax.set_zticks([])
ax.set_title("Branch from Neuron Trace")
plt.axis('on')
plt.show()
../../_images/notebooks_algorithms_biccn_demo_10_0.png
[6]:
curvature = spline_fxns.curvature(u, tck[0], tck[1], tck[2])
ra_eq, cdf = hist_equalize(curvature, 100)

fig = plt.figure(figsize=(12, 10), dpi=80)
ax = fig.add_subplot(111,projection="3d")
im = ax.scatter(pts[0], pts[1], pts[2], c=curvature, cmap='Reds')
ax.set_xlabel('x ($\mu m$)')
ax.set_ylabel('y ($\mu m$)')
ax.set_zlabel('z ($\mu m$)')
ax.set_title("Branch from Neuron Trace Showing Curvature")
cbar = fig.colorbar(im)
cbar.set_label('Curvature ($(\mu m)^{-1}$)')
../../_images/notebooks_algorithms_biccn_demo_11_0.png
[7]:
ra_eq, cdf = hist_equalize(curvature, 100)

fig = plt.figure(figsize=(12, 10), dpi=80)
ax = fig.add_subplot(111,projection="3d")
im = ax.scatter(pts[0], pts[1], pts[2], c=ra_eq, cmap='Reds')
ax.set_xlabel('x ($\mu m$)')
ax.set_ylabel('y ($\mu m$)')
ax.set_zlabel('z ($\mu m$)')
ax.set_title("Branch from Neuron Trace Showing Curvature Percentile")
cbar = fig.colorbar(im)
cbar.set_label('Curvature Percentile')
../../_images/notebooks_algorithms_biccn_demo_12_0.png
[8]:
torsion = spline_fxns.torsion(u, tck[0], tck[1], tck[2])
ra_eq, cdf = hist_equalize(torsion, 100)

fig = plt.figure(figsize=(12, 10), dpi=80)
ax = fig.add_subplot(111,projection="3d")
im = ax.scatter(pts[0], pts[1], pts[2], c=ra_eq, cmap='Reds')
ax.set_xlabel('x ($\mu m$)')
ax.set_ylabel('y ($\mu m$)')
ax.set_zlabel('z ($\mu m$)')
ax.set_title("Branch from Neuron Trace Showing Torsion Percentile")
cbar = fig.colorbar(im)
cbar.set_label('Torsion Percentile')
../../_images/notebooks_algorithms_biccn_demo_13_0.png
[ ]: