In [108]:
import sys
sys.path.append('../code/functions/')
import tiffIO as io
import math
import cv2
import time
import pickle
import numpy as np
import synapseLib as sl
import matplotlib.pyplot as plt
from skimage.exposure import equalize_adapthist
from scipy.ndimage.filters import convolve
from skimage.filters import sobel
from skimage.morphology import dilation
from scipy.spatial import KDTree
from random import randint
from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot
import plotly.graph_objs as go
init_notebook_mode(connected=True)
In [2]:
data = np.array(io.loadTiff('../data/rr46b_s0_ch1.tif'))
In [3]:
axons, _ = sl.extractAxons(np.stack([equalize_adapthist(elem) for elem in data]), percentile=50)
In [18]:
plt.imshow(axons[15], cmap='gray')
plt.show()
In [76]:
img = data[15]
for kernelSize in [8]:
img = convolve(img, np.ones((kernelSize, kernelSize)))
img = equalize_adapthist(img)
plt.imshow(img, cmap='gray')
plt.show()
In [61]:
axons2, _ = sl.extractAxons([img], percentile=50)
In [63]:
plt.imshow(axons2[0], cmap='gray')
plt.show()
In [4]:
def evolveAxons(data, epochs=2, n=2, neighborhood=16, dilations=5, percentile=50):
species = data
kernel = np.ones((8, 8))
for i in range(epochs):
print i
genus, _ = sl.extractAxons(np.stack([equalize_adapthist(elem) for elem in species]),
percentile=percentile,
neighborhood=neighborhood,
n = n,
dilations=dilations)
species = np.stack([convolve(elem, kernel) for elem in genus]).astype(np.int64)
return species
In [5]:
axons = evolveAxons(data[10:20], epochs=10, dilations=10)
In [138]:
plt.imshow(axons[5], cmap='gray')
plt.show()
In [141]:
kernelX = [[1,0,-1],
[2, 0, -2],
[1, 0, -1]]
kernelY = [[1, 2, 1],
[0, 0, 0],
[-1, -2, -1]]
In [143]:
xGrad = convolve(axons[5], kernelX)
yGrad = convolve(axons[5], kernelY)
grad = np.sqrt(np.add(np.power(xGrad, 2), np.power(yGrad, 2)))
In [144]:
plt.imshow(grad, cmap='gray')
plt.show()
In [161]:
kernel = np.ones((16, 16))
kernel[1:15, 1:15] = 0
'''
kernel[7, 7] = -6
kernel[7, 8] = -6
kernel [8, 7] = -6
kernel[8, 8] = -6
'''
Out[161]:
In [162]:
nodes = convolve(grad, kernel)
plt.imshow(nodes, cmap='gray')
plt.show()
In [166]:
test = np.logical_and(np.logical_xor(nodes, grad), axons[5])
plt.imshow(test, cmap='gray')
plt.show()
In [8]:
def generateNodeImg(axonImg, step=64):
kernelX = [[1,0,-1],
[2, 0, -2],
[1, 0, -1]]
kernelY = [[1, 2, 1],
[0, 0, 0],
[-1, -2, -1]]
xGrad = convolve(axons[5], kernelX)
yGrad = convolve(axons[5], kernelY)
grad = np.sqrt(np.add(np.power(xGrad, 2), np.power(yGrad, 2)))
symmetryKernel = np.ones((16, 16))
symmetryKernel[1:15, 1:15] = 0
symmetryKernel[7:9, 7:9] = -1
potentialNodes = convolve(grad, symmetryKernel)
nodes = np.multiply(np.logical_and(np.logical_xor(potentialNodes, grad), axonImg), axonImg)
nonMaxSuppression = np.zeros_like(nodes)
for y in range(0, 1024, step):
for x in range(0, 1024, step):
sub = nodes[y:y+step, x:x+step]
aMax = np.argmax(sub)
yMax = aMax/step
xMax = aMax%step
nonMaxSuppression[y+yMax, x+xMax] = sub[yMax, xMax]
return nonMaxSuppression
In [94]:
nodeImg = generateNodeImg(axons[5], 64)
plt.imshow(nodeImg, cmap='gray')
plt.show()
In [95]:
plt.imshow(axons[5], cmap='gray')
plt.show()
In [96]:
print np.count_nonzero(nodeImg)
print len(zip(*(np.nonzero(nodeImg))))
In [97]:
axCp = axons[5].copy()
for node in zip(*(np.nonzero(nodeImg))):
axCp[node[0]-5:node[0]+5, node[1]-5:node[1]+5] = 255
plt.imshow(axCp, cmap='gray')
plt.show()
In [104]:
def generateGraph(nodeImg, axons, thickness=10, meanThresh = 50, devThresh = 10):
axCp = axons.copy()
aves = []
devs = []
potEdges = []
edges = []
nodes = zip(*(np.nonzero(nodeImg)))
for i in range(len(nodes)):
print i/float(len(nodes))
for j in range(i+1, len(nodes)):
y0, x0 = nodes[i]
y1, x1 = nodes[j]
length = int(np.hypot(x1-x0, y1-y0))
x, y = np.linspace(x0, x1, length).astype(int), np.linspace(y0, y1, length).astype(int)
potEdges.append([y, x])
potEdgeStats = []
for k in range(length):
sub = axons[max(y[k]-thickness,0):min(y[k]+thickness, 1024), max(x[k]-thickness, 0):min(x[k]+thickness, 1024)]
potEdgeStats.append(np.mean(sub))
aves.append(np.mean(potEdgeStats))
devs.append(np.std(potEdgeStats))
meanCut = np.percentile(aves, meanThresh)
devCut = np.percentile(devs, devThresh)
for i in range(len(potEdges)):
if aves[i] > meanCut and devs[i] < devCut:
edges.append(potEdges[i])
for edge in edges:
y = edge[0]
x = edge[1]
for k in range(len(y)):
axCp[max(y[k]-3,0):min(y[k]+3, 1024), max(x[k]-3, 0):min(x[k]+3, 1024)] = 255
return nodes, edges, axCp
In [105]:
nodes, edges, vis = generateGraph(nodeImg, axons[5])
In [106]:
plt.imshow(vis, cmap='gray')
plt.show()
In [157]:
def estimateGraph(nodeImg, axons, thickness=10, neighbors=6, baselineSize=10):
axCp = axons.copy()
baseline = []
for i in range(baselineSize):
y0, x0 = randint(0, axons.shape[0]), randint(0, axons.shape[1])
y1, x1 = randint(0, axons.shape[0]), randint(0, axons.shape[1])
length = int(np.hypot(x1-x0, y1-y0))
x, y = np.linspace(x0, x1, length).astype(int), np.linspace(y0, y1, length).astype(int)
edgeStats = []
for k in range(length):
sub = axons[max(y[k]-thickness,0):min(y[k]+thickness, 1024), max(x[k]-thickness, 0):min(x[k]+thickness, 1024)]
edgeStats.append(np.mean(sub))
baseline.append(np.mean(edgeStats))
baseMu = np.mean(baseline)
baseSig = np.std(baseline)
edges = []
nodes = zip(*(np.nonzero(nodeImg)))
tree = KDTree(nodes)
for curIdx, node in enumerate(nodes):
print curIdx/float(len(nodes))
partnerIdxList = tree.query(node, k=neighbors)[1]
partners = []
for partnerIdx in partnerIdxList:
if partnerIdx > curIdx:
partners.append(nodes[partnerIdx])
for partner in partners:
y0, x0 = node
y1, x1 = partner
length = int(np.hypot(x1-x0, y1-y0))
x, y = np.linspace(x0, x1, length).astype(int), np.linspace(y0, y1, length).astype(int)
edgeStats = []
for k in range(length):
sub = axons[max(y[k]-thickness,0):min(y[k]+thickness, 1024), max(x[k]-thickness, 0):min(x[k]+thickness, 1024)]
edgeStats.append(np.mean(sub))
dp = np.mean(edgeStats)
z = (dp - baseMu)/float(baseSig)
if z > 1.5:
edges.append([y, x])
for edge in edges:
y = edge[0]
x = edge[1]
for k in range(len(y)):
axCp[max(y[k]-3,0):min(y[k]+3, 1024), max(x[k]-3, 0):min(x[k]+3, 1024)] = 255
return nodes, edges, axCp
In [158]:
nodes, edges, vis = estimateGraph(nodeImg, axons[5])
In [160]:
class anatomy:
def __init__(self, nodes, edges, vis):
self._nodes = nodes
self._edges = edges
self._vis = vis
def generateAnatomyVolume(data):
anatomyVolume = []
axonImgs = evolveAxons(data, epochs=10, dilations=10)
for axonImg in axonImgs:
nodeImg = generateNodeImg(axonImg, 64)
nodes, edges, vis = estimateGraph(nodeImg, axonImg)
anatomyVolume.append(anatomy(nodes, edges, vis))
return anatomyVolume
In [167]:
data2 = np.array(io.loadTiff('../data/rr46b_s1_ch1.tif'))
start = time.time()
anatomyVol = generateAnatomyVolume(data2[10:15])
In [181]:
plt.imshow(equalize_adapthist(data2[10])*1000, cmap='gray')
plt.show()
In [174]:
plt.imshow(anatomyVol[0]._vis, cmap='gray')
plt.show()
In [159]:
plt.imshow(vis, cmap='gray')
plt.show()
In [ ]: