In [5]:
import sys
sys.path.append('../code/functions/')
import tiffIO as io
import neuroGraphLib as ngl
from neuroGraphLib import neuroGraph
from scipy.stats import zscore
from skimage.exposure import equalize_adapthist
from scipy.ndimage import convolve
from scipy.signal import argrelmax
import synapseLib as sl
from skimage.measure import label
import matplotlib.pyplot as plt
import cv2
import time
import math
import numpy as np
In [6]:
data_0 = np.array(io.loadTiff('../data/rr46b_s0_ch1.tif'))
data_1 = np.array(io.loadTiff('../data/rr46b_s1_ch1.tif'))
In [3]:
start = time.time()
neuroGraphStack_0 = ngl.generateNeuroGraphStack(data_0)
end = time.time()
print 'Neuro Graph Generation Took: ', end - start
In [7]:
start = time.time()
neuroGraphStack_1 = ngl.generateNeuroGraphStack(data_1)
end = time.time()
print 'Neuro Graph Generation Took: ', end - start
In [5]:
for i in range(10, 20):
plt.subplot(121)
plt.imshow(neuroGraphStack_0[i]._vis, cmap='gray')
plt.subplot(122)
plt.imshow(data_0[i]*100, cmap='gray')
plt.show()
In [6]:
for i in range(10, 20):
plt.subplot(121)
plt.imshow(neuroGraphStack_1[i]._vis, cmap='gray')
plt.subplot(122)
plt.imshow(data_1[i]*100, cmap='gray')
plt.show()
In [55]:
def generateTube(mat, x0, x1, y0, y1, thickness):
x = np.linspace(x0, x1, int(np.sqrt(np.abs(x1**2 - x0**2))))
y = np.linspace(y0, y1, int(np.sqrt(np.abs(y1**2 - y0**2))))
coords = zip(x, y)
perpSlope = -1.0*(x1 - x0)/(y1 - y0)
numSamplePoints = thickness
netCoords = []
for coord in coords:
xtop = int(math.floor(coord[0])) + numSamplePoints*perpSlope
xbot = int(math.floor(coord[0])) - numSamplePoints*perpSlope
xperp = np.linspace(xtop, xbot, int(np.sqrt(np.abs(xtop**2 - xbot**2))))
ytop = int(math.floor(coord[1])) + numSamplePoints
ybot = int(math.floor(coord[1])) - numSamplePoints
yperp = np.linspace(ytop, ybot, int(np.sqrt(np.abs(ytop**2 - ybot**2))))
coordsPerp = zip(xperp, yperp)
finalCoordsPerp = []
for coordPerp in coordsPerp:
if coordPerp[0] > 0 and coordPerp[1] > 0 and coordPerp[0] < mat.shape[0] - 1 and coordPerp[1] < mat.shape[1] - 1:
finalCoordsPerp.append(((int(math.floor(coordPerp[1]))), (int(math.floor(coordPerp[0])))))
netCoords.append(finalCoordsPerp)
return netCoords
mat = np.zeros((301,301))
netCoords = generateTube(mat, 0, 300, 0, 300, 50)
#finalCoordsPerp is what Brandon wants
for viewCoord in netCoords:
for subCoord in viewCoord:
mat[subCoord[0]][subCoord[1]] = 2
plt.imshow(mat)
plt.show()
In [8]:
testNeuroGraph = neuroGraphStack_1[15]
In [57]:
len(testNeuroGraph._edges)
Out[57]:
In [58]:
vis = np.zeros_like(testNeuroGraph._vis)
for i in range(0, len(testNeuroGraph._edges), 10):
y, x = testNeuroGraph._edges[i]
curROI = generateTube(vis, x[0], x[-1], y[0], y[-1], 10)
ROIDist = []
for orthagonal in curROI:
orthagonalDist = []
for elem in orthagonal:
orthagonalDist.append(data_1[15][elem[0]][elem[1]])
ROIDist.append(np.mean(orthagonalDist))
for k in range(len(y)):
vis[max(y[k]-3,0):min(y[k]+3, vis.shape[0]), max(x[k]-3, 0):min(x[k]+3, vis.shape[1])] = 255
plt.figure()
plt.scatter(range(len(ROIDist)), ROIDist)
plt.show()
In [ ]:
plt.figure()
plt.imshow(vis, cmap='gray')
plt.show()
As expected, the distribution spikes when walking across the edge. Time to mark the elements that are sufficiently 'spiked' as synapse
In [24]:
def getSynapseROIs(neuroGraph):
synapseROIs = []
for i in range(0, len(testNeuroGraph._edges)):
print i/float(len(testNeuroGraph._edges))
y, x = testNeuroGraph._edges[i]
start = time.time()
curROI = generateTube(vis, x[0], y[0], x[-1], y[-1], 10)
end = time.time()
print 'Generating tube took: ', end - start
start = time.time()
ROIDist = []
for orthagonal in curROI:
orthagonalDist = []
for elem in orthagonal:
orthagonalDist.append(data_1[20][elem[0]][elem[1]])
ROIDist.append(np.mean(orthagonalDist))
end = time.time()
print 'Generating dist took: ', end-start
start = time.time()
mu = np.average(ROIDist)
sig = np.std(ROIDist)
end = time.time()
print'Numpy stats took: ', end - start
start = time.time()
z = [(elem - mu)/sig for elem in ROIDist]
for i, elem in enumerate(z):
if z >= 2:
synapseROIs = synapseROIs + curROI[i]
end = time.time()
print 'List stats took: ', end - start
return synapseROIs
In [25]:
synapseROIs = getSynapseROIs(testNeuroGraph)
In [59]:
def getSynapseROIs(neuroGraph):
synapseROIs = []
#for i in range(0, 1):
for i in range(0, len(testNeuroGraph._edges)):
print i/float(len(testNeuroGraph._edges))
try:
y, x = testNeuroGraph._edges[i]
curROI = np.array(generateTube(vis, x[0], x[-1], y[0], y[-1], 20))
ROIDist = []
for orthagonal in curROI:
orthagonalDist = []
for elem in orthagonal:
#since will returns xy
orthagonalDist.append(data_1[20][elem[0]][elem[1]])
ROIDist.append(np.mean(orthagonalDist))
valid = np.array(np.where(zscore(ROIDist) > 1.5)[0])
if len(valid) > 0:
toAdd = np.concatenate([elem for elem in curROI[valid]])
synapseROIs.append(toAdd)
except:
continue
return synapseROIs
In [60]:
synapseROIs = getSynapseROIs(testNeuroGraph)
In [ ]:
roi = testNeuroGraph._vis.copy()
for orth in synapseROIs:
for y, x in orth:
roi[y][x] = 255
plt.imshow(roi, cmap='gray')
plt.show()
In [9]:
def generateTube(x0, x1, y0, y1, thickness, mat):
x = np.linspace(x0, x1, int(np.sqrt(np.abs((x1 - x0)**2 + (y1 - y0)**2))))
y = np.linspace(y0, y1, int(np.sqrt(np.abs((x1 - x0)**2 + (y1 - y0)**2))))
coords = zip(x, y)
perpSlope = None
if y1-y0 == 0:
perpSlope = 1000
elif x1-x0 == 0:
perpSlope == .0001
else:
perpSlope = -1.0*(x1 - x0)/(y1 - y0)
numSamplePoints = thickness/2
netCoords = []
for coord in coords:
xtop = coord[0] + numSamplePoints/np.sqrt(1 + perpSlope*perpSlope)
xbot = coord[0] - numSamplePoints/np.sqrt(1 + perpSlope*perpSlope)
ytop = coord[1] + perpSlope*numSamplePoints/np.sqrt(1 + perpSlope*perpSlope)
ybot = coord[1] - perpSlope*numSamplePoints/np.sqrt(1 + perpSlope*perpSlope)
xperp = np.linspace(xtop, xbot, int(np.sqrt(np.abs((xtop - xbot)**2 + (ytop - ybot)**2))))
yperp = np.linspace(ytop, ybot, int(np.sqrt(np.abs((xtop - xbot)**2 + (ytop - ybot)**2))))
coordsPerp = zip(xperp, yperp)
finalCoordsPerp = []
for coordPerp in coordsPerp:
if coordPerp[0] > 0 and coordPerp[1] > 0 and coordPerp[0] < mat.shape[0] - 1 and coordPerp[1] < mat.shape[1] - 1:
finalCoordsPerp.append(((int(math.floor(coordPerp[1]))), (int(math.floor(coordPerp[0])))))
netCoords.append(finalCoordsPerp)
return netCoords
In [10]:
def waveletSynapseROIs(neuroGraph):
boostedImg, _ = ngl.boostDendrites(data_1[13:17], percentile=50)
adaptiveImg = ngl.adaptiveThreshold(boostedImg, 32, 32)[2]
plt.imshow(adaptiveImg, cmap='gray')
plt.show()
synapseROIs = []
#for i in range(0, 1):
for i in range(0, len(testNeuroGraph._edges)):
print i/float(len(testNeuroGraph._edges))
try:
y, x = testNeuroGraph._edges[i]
curROI = np.array(generateTube(x[0], x[-1], y[0], y[-1], 30, adaptiveImg))
ROIDist = []
ROIDistLedger = []
for orthagonal in curROI:
orthagonalDist = []
orthagonalMaxLedger = []
for elem in orthagonal:
#since will returns xy
orthagonalDist.append(adaptiveImg[elem[0]][elem[1]])
orthagonalMaxLedger.append((elem[0], elem[1]))
ROIDist.append(np.max(orthagonalDist))
ROIDistLedger.append(orthagonalMaxLedger[np.argmax(orthagonalDist)])
valid = np.array(np.where(zscore(ROIDist) > 1)[0])
'''
valid = argrelmax(np.array(ROIDist))[0]
if len(valid) > 0:
toAdd = np.concatenate([elem for elem in curROI[valid]])
synapseROIs.append(toAdd)
if len(valid) > 0:
toAdd = np.concatenate([elem for elem in curROI[valid]])
synapseROIs.append(toAdd)
'''
if len(valid) > 0:
for idx in valid:
synapseROIs.append(ROIDistLedger[idx])
else:
print 'no maxima found!'
except:
print sys.exc_info()[0]
return synapseROIs
In [11]:
synapseROIs = waveletSynapseROIs(testNeuroGraph)
In [12]:
roi = testNeuroGraph._vis.copy()
#for orth in synapseROIs:
for i in range(1):
for y, x in synapseROIs: # was orth
roi[y][x] = 9001
plt.imshow(roi, cmap='gray')
plt.show()
In [17]:
data_1_b, _ = ngl.boostDendrites(data_1)
print data_1_b.shape
zoom = data_1_b[15, 600:1000, 300:700]
plt.imshow(zoom, cmap='gray')
plt.show()
roi = np.zeros_like(data_1[15].copy())
#for orth in synapseROIs:
for i in range(1):
for y, x in synapseROIs: # was orth
roi[y][x] = 300
kernel_a = np.ones((8, 8))
kernel_b = np.ones((3, 3))
roi_c = np.add(convolve(roi, kernel_a), convolve(roi, kernel_b))
plt.imshow(roi_c[100:400, 500:900], cmap='gray')
#plt.imshow(roi, cmap='gray')
plt.show()
In [69]:
dat = np.extract(roi_c[600:1000, 300:700] > 0, roi_c)
plt.hist(dat, bins=50)
plt.show()
In [ ]: